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£N| ' Abstract 

These are the notes for a set of lectures delivered by the two authors at the Les Houches Summer 
School on 'Complex Systems' in July 2006. They provide an introduction to the basic concepts in modern 
(probabilistic) coding theory, highlighting connections with statistical mechanics. We also stress common 
concepts with other disciplines dealing with similar problems that can be generically referred to as 'large 
graphical models'. 

While most of the lectures are devoted to the classical channel coding problem over simple memoryless 
channels, we present a discussion of more complex channel models. We conclude with an overview of 

■ the main open challenges in the field. 

m 

. 1 Introduction and Outline 

The last few years have witnessed an impressive convergence of interests between disciplines which are 
a priori well separated: coding and information theory, statistical inference, statistical mechanics (in 
particular, mean field disordered systems), as well as theoretical computer science. The underlying 
reason for this convergence is the importance of probabilistic models and/or probabilistic techniques in 
each of these domains. This has long been obvious in information theory [53], statistical mechanics [TP] , 
and statistical inference [IS]. In the last few years it has also become apparent in coding theory and 
theoretical computer science. In the first case, the invention of Turbo codes [7] and the re-invention 
of Low-Density Parity-Check (LDPC) codes [30, [28] has motivated the use of random constructions for 
coding information in robust/compact ways [50]. In the second case (theoretical computer science) the 
relevance of randomized algorithms has steadily increased (see for instance [41]). thus motivating deep 
theoretical developments. A particularly important example is provided by the Monte Carlo Markov 
Chain method for counting and sampling random structures. 

Given this common probabilistic background, some analogies between these disciplines is not very 
surprising nor is it particularly interesting. The key new ingredient which lifts the connections beyond 
some superficial commonalities is that one can name specific problems, questions, and results which lie 
at the intersection of these fields while being of central interest for each of them. The set of problems 
and techniques thus defined can be somewhat loosely named "theory of large graphical models." The 
typical setting is the following: a large set of random variables taking values in a finite (typically quite 
small) alphabet with a "local" dependency structure; this local dependency structure is conveniently 
described by an appropriate graph. 



*The work of A. Montanari was partially supported by the European Union under the project EVERGROW. The work of 
R. Urbanke was partially supported by the NCCR-MICS, a center supported by the Swiss National Science Foundation under 
grant number 5005-67322. 
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In this lecture we shall use "modern" coding theory as an entry point to the domain. There are several 
motivations for this: (i) theoretical work on this topic is strongly motivated by concrete and well-defined 
practical applications; (ii) the probabilistic approach mentioned above has been quite successful and has 
substantially changed the field (whence the reference to modern coding theory); (in) a sufficiently 
detailed picture exists illustrating the interplay among different view points. 

We start in Section[2]with a brief outline of the (channel coding) problem. This allows us to introduce 
the standard definitions and terminology used in this field. In Section [3] we introduce ensembles of 
codes defined by sparse random graphs and discuss their most basic property - the weight distribution. 
In Section [4] we phrase the decoding problem as an inference problem on a graph and consider the 
performance of the efficient (albeit in general suboptimal) message-passing decoder. We show how the 
performance of such a combination (sparse graph code and message-passing decoding) can be analyzed 
and we discuss the relationship of the performance under message-passing decoding to the performance 
of the optimal decoder. In Section [5] we briefly touch on some problems beyond coding, in order to 
show as similar concept emerge there. In particular, we discuss how message passing techniques can be 
successfully used in some families of counting/inference problems. In Section [6] we show that several of 
the simplifying assumptions (binary case, symmetry of channel, memoryless channels) are convenient in 
that they allow for a simple theory but are not really necessary. In particular, we discuss a simple channel 
with memory and we see how to proceed in the asymmetric case. Finally, we conclude in Section [7] with 
a few fundamental open problems. 

To readers who would like to find current contributions on this topic we recommend the IEEE 
Transactions on Information Theory. A considerably more in-depth discussion can be found in the two 
upcoming books Information, Physics and Computation [36j and Modern Coding Theory |50j . Standard 
references on coding theory are [6] [9] [26] and very readable introductions to information theory can be 
found in [121 ED] . Other useful reference sources are the book by Nishimori [44] as well as the book by 
MacKay [29]. 

2 Background: The Channel Coding Problem 

The central problem of communications is how to transmit information reliably through a noisy (and 
thus unreliable) communication channel. Coding theory aims at accomplishing this task by adding a 
properly designed redundancy to the transmitted message. This redundancy is then used at the receiver 
to reconstruct the original message despite the noise introduced by the channel. 

2.1 The Problem 

In order to model the situation described above we shall assume that the noise is random with some 
known distribution^] To keep things simple we shall assume that the communication channel admits as 
input binary symbols x € {0, 1}, while the output belongs to some finite alphabet A. We denote the 
probability of observing the output y £ A given that the input was x G {0, 1} by Q(y\x). The channel 
model is defined by the transition probability matrix 

Q = {Q{y\x) : x e {0, 1}, y e A} . (2.1) 

Of course, the entries of this matrix must be non-negative and normalized in such a way that ^ Q(y\x) = 
1. It is convenient to have a few simple examples in mind. We refer to Fig. [T]for an illustration of the 
channel models which we introduce in the following three examples. 

Example 1: The binary symmetric channel BSC(p) is defined by letting A — {0, 1} and Q(0|0) = 
Q(l|l) = 1 — p; the normalization then enforces Q(1|0) = Q(0| 1) = p. In words, the channel "flips" 
the input bit with probability p G [0, 1]. Since flips are introduced for each bit independently we 



1 lt is worth mentioning that an alternative approach would be to consider the noise as 'adversarial' (or worst case) under 
some constraint on its intensity. 
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Figure 1: Schematic description of three simple binary memory less channels. From left to right: binary 
symmetric channel BSC(p), binary erasure channel BEC(e), and Z channel ZC(p). 



say that the channel is memoryless. Except for an example in Section 16.21 all channels which we 
consider are memoryless. 

Example 2: The binary erasure channel BEC(e) is defined by A = {0, 1, *} and Q(0|0) = Q(l|l) = 
1 — e while Q(*\0) — Q(*|l) = e. In words, the channel input is erased with probability e and it is 
transmitted correctly otherwise. 

Example 3: The Z-channel ZC(p) has an output alphabet A = {0, 1} but acts differently on input 
(that is transmitted correctly) and 1 (that is flipped with probability p) . We invite the reader 
to write the transition probability matrix. 

Since in each case the input is binary we speak of a binary-input channel. Since further in all models each 
input symbol is distorted independently from all other ones we say that the channels are memoryless. 
It is convenient to further restrict our attention to symmetric channels: this means that there is an 
involution on A (i.e. a mapping t ; A — ► A such that ioi = 1) so that Q(y\0) = Q(b{y)\l). (E.g., 
if A = K then we could require that Q(y|0) = Q(—y\l).) This condition is satisfied by the first two 
examples above but not by the third one. To summarize these three properties one refers to such models 
as BMS channels. 

In order to complete the problem description we need to formalize the information which is to 
be transmitted. We shall model this probabilistically as well and assume that the transmitter has 
an information source that provides an infinite stream of i.i.d. fair coins: {zf,i — 0,1,2,...}, with 
Zi G {0, 1} uniformly at random. The goal is to reproduce this stream faithfully after communicating it 
over the noisy channel. 

Let us stress that, despite its simplification, the present setting contains most of the crucial and 
challenges of the channel coding problem. Some of the many generalizations are described in Section [6l 

2.2 Block Coding 

The (general) coding strategy we shall consider here is block coding. It works as follows: 

• The source stream {zi} is chopped into blocks of length L. Denote one such block by z, z = 
{z u ...,z L ) G {0,1} L . 

• Each block is fed into an encoder. This is a map F : {0, 1} L — > {0, 1}^, for some fixed N > L (the 
blocklength) . In words, the encoder introduces redundancy in the source message. Without loss of 
generality we can assume F to be injective. It this was not the case, even in the absence of noise, 
we could not uniquely recover the transmitted information from the observed codeword. 

• The image of {0, 1} L under the map F is called the codebook, or sometimes the code, and it will 
be denoted by £. The code contains \€\ = 2 L strings of length N called codewords. These are the 
possible channel inputs. The codeword x = F(z) is sent through the channel, bit by bit. 
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Figure 2: Flow chart of a block coding scheme. 



• Let y — (yi, ■ • • , j/jv) € -A be the channel output. Conditioned on x the t/,-, i = 1, • • • , L, are 
independent random variables with distribution yi = Q(- \xi) (here and below = denotes identity 
in distribution and x = P( • ) means that a; is a random variable with distribution P( • )). 

• The channel output is fed into a decoder, which is a map F : A — ► {0, 1} L . It is the objective of 
the decoder to reconstruct the source z from the noisy channel output y. 

The flow chart describing this coding scheme is shown in Fig. O It is convenient to slightly modify 
the above scheme. Notice that, under the hypothesis that the encoder is injective, the codebook is in 
one-to-one correspondence with the source sequences. Since these are equiprobable, the transmitted 
codewords are equiprobable as well. We can therefore equivalently assume that the transmitter picks a 
codeword uniformly at random and transmits it. Every reference to the source stream can be eliminated 
if we redefine the decoder to be a map F : A N — > {0, 1}^, i.e., the decoder aims to reconstruct the 
transmitted codeword. If F(y) £ £ we declare an error @ In the following we shall also use the notation 

F(y) =Mv) = (xi(y),---,x N (y)). 

One crucial parameter of a code is its rate: it quantifies how many bits of information are transmitted 
per channel use, 

i^£ = ilog 2 |£|. (2.2) 

Two fundamental performance parameters are the bit (or 'symbol') and block (or 'word') error rates. 
The block error rate is the probability that the input codeword is not recovered correctly at the end of 
the process, 

Pb = P {My) ? x) . (2.3) 
The bit error rate is the expected fraction of bits that are not recovered correctly, 

N 

X 



1 N 

Pb^E P {^)^} ■ ( 2 - 4 ) 



It should not be too surprising that one can trade-off rate and error probability. We want to achieve 
a high rate and achieve a low probability of error. However, increasing the rate decreases the redundancy 
built into the codeword, thus inducing a higher error probability. The aim of coding theory is to choose 
the code € and the decoding function x( • ) in a way to optimize this trade-off. 



2 More precisely, if we are interested only in the block probability of error, i.e., the frequency at which the whole block of 
data is decoded correctly, then indeed any one-to-one mapping between information word and codeword performs identical. If, 
on the other hand, we are interested in the fraction of bits that we decode correctly then the exact mapping from information 
word to codeword does come into play. We shall ignore this somewhat subtle point in the sequel. 
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2.3 Decoding 

Given the code there is a simple (although in general not computationally efficient) prescription for the 
decoder. If we want to minimize the block error rate, we must chose the most likely codeword, 

x B (y) = argmaxP{A = x\Y = y) . (2.5) 

To minimize the bit error rate we must instead return the sequence of most likely bits, 

x\ (y) = arg max F{X t = Xi \Y = y) . (2.6) 

The reason of these prescriptions is the object of the next exercise. 

Exercise 1: Let (U, V) be a pair of discrete random variables. Think of U as a 'hidden' variable 
and imagine you observe V = v. We want to understand what is the optimal estimate for U given 
V — v. Show that the function v <— > u(v) that minimizes the error probability P(u) = P {U ^ u(V)} 
is given by 

u(v) = argmax F {U = u\V = v} . (2.7) 

u 

It is instructive to explicitly write down the conditional distribution of the channel input given the 
output. We shall denote it as fJ.<r, y (x) = P{X_ = x\Y_ = y} (and sometimes drop the subscripts £ and y 
if they are clear from the context). Using Bayes rule we get 

1 N 

Ht,v(x) = Y(£y) n^l^Iefe), (2-8) 

where denotes the code membership function (Ic(x) = 1 if x € £ and = otherwise). 

According to the above discussion, decoding amounts to computing the marginals (for symbol MAP) 
or the model (for word MAP) of n(-). More generally, we would like to understand the properties of 
/x( • ): is it concentrated on a single codeword or spread over many of them? In the latter case, are these 
close to each other or very different? And what is their relationship with the transmitted codeword? 

The connection to statistical mechanics emerges in the study of the decoding problem 15"T] . To 
make it completely transparent we rewrite the distribution /i( • ) in Boltzmann form 

E , x) = ( -E^logQte), ifaec, (21Q) 

[ +oo, otherwise. 
The word MAP and bit MAP rule can then be written as 

x B (y) = argmin£ , £,„(x) , (2-11) 

— X 

%i{y) = argmax V (i<r, y (x) . (2.12) 

— Xi — ' 

In words, word MAP amounts to computing the ground state of a certain energy function, and bit 
MAP corresponds to computing the expectation with respect to the Boltzmann distribution. Notice 
furthermore that /i( • ) is itself random because of the randomness in y (and we shall introduce further 
randomness in the choice of the code). This is analogous to what happens in statistical physics of 
disordered systems, with y playing the role of quenched random variables. 



We recall that the mode of a distribution with density /i( • ) is the value of x that maximizes fi(x). 
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2.4 Conditional Entropy and Free Energy 

As mentioned above, we are interested in understanding the properties of the (random) distribution 
/•*£,»( ' )■ O ne possible way of formalizing this idea is to consider the entropy of this distribution. 

Let us recall that the (Shannon) entropy of a discrete random variable X (or, equivalently, of its 
distribution) quantifies, in a very precise sense, the 'uncertainty' associated with A0 It is given by 

H(X) = -J2v{x)logF(x) . (2.13) 

X 

For two random variables X and Y one defines the conditional entropy of X given Y as 

H{X\Y) = -^¥{x,y)\ogf>{x\y)=¥< y \-^Y(x\Y)\ogI>{x\Y)\. (2.14) 

x,y \ x ) 

This quantifies the remaining uncertainty about X when Y is observed. 

Considering now the coding problem. Denote by X_ the (uniformly random) transmitted codeword 
and by Y_ the channel output. The right-most expression in Eq. (|2.14p states that H(X\Y) is the 
expectation of the entropy of the conditional distribution ^c,y{ ' ) with respect to y. 

Let us denote by v<t{x) the probability that a uniformly random codeword in £ takes the value x at 
the i-th position, averaged over i. Then a straightforward calculation yields 

1 N { 1 N 1 

H(x\Y) = -—^rjg^i^iogj^-^n^l^)] > ( 2 - 15 ) 

= -Nj2Mx)Q(y\x)logQ(y\x)+EylogZ(€,y). (2.16) 

x,y 

The 'type' u<t{x) is usually a fairly straightforward characteristic of the code. For most of the examples 
considered below we can take v^{0) — v<t(l) = 1/2. As a consequence the first of the terms above is 
trivial to compute (it requires summing over 2\A\ terms). 

On the other hand the second term is highly non-trivial. The reader will recognize the expectation 
of a free energy, with y playing the role of a quenched random variable. 

The conditional entropy H(X_\Y_) provides an answer to the question: how many codewords is /J-c. y ( ■ ) 
spread over? It turns out that about e H (— '— ' of them carry most of the weight. 

2.5 Shannon Theorem and Random Coding 

As mentioned above, there exists an obvious tradeoff between high rate and low error probability. In his 
celebrated 1948 paper [53], Shannon derived the optimal error probability- vs-rate curve in the limit of 
large blocklengths. In particular, he proved that if the rate is larger than a particular threshold, then 
the error probability can be made arbitrarily small. The threshold depends on the channel and it is 
called the channel capacity. The capacity of a BMS channel (measured in bits per channel use) is given 
by the following elementary expression, 

C(Q) =H(X)-H(X\Y) 

For instance, the capacity of a BSC(p) is C(p) = 1 — f)2(p), (where f)2(p) = — plog 2 p— (1— p) log 2 (l — p) is 
the entropy of a Bernoulli random variable of parameter p) while the capacity of a BEC(e) is C(e) = 1 — e. 
As an illustration, the capacity of a BSC(p) with flip probability p ~ 0.110028 is C(p) = 1/2: such a 
channel can be used to transmit reliably 1/2 bit of information per channel use. 



4 For a very readable account of information theory we recommend [12] , 
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Theorem 2.1 (Channel Coding Theorem). For any BMS channel with transition probability Q and 
R < C(Q) there exists a sequence of codes Cat of increasing blocklength N and rate Rn — * R whose block 
error probability Pg^ — > as N — > oo. 

Fice versa, for any R > C(Q) the block error probability of a code with rate at least R is bounded 
away from 0. 

The prove of the first part ('achievability') is one of the first examples of the so-called 'probabilistic 
method'. In order to prove that there exists an object with a certain property (a code with small 
error probability), one constructs a probability distribution over all potential candidates (all codes 
of a certain blocklength and rate) and shows that a random element has the desired property with 
non-vanishing probability. The power of this approach is in the (meta-mathematical) observation that 
random constructions are often much easier to produce than explicit, deterministic ones. 
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Figure 3: Exponential growth rate for the expected distance enumerator KN x (o)(n5) within the random 
code ensemble. 



The distribution over codes proposed by Shannon is usually referred to as the random code (or, 
Shannon) ensemble, and is particularly simple. One picks a code £ uniformly at random among all codes 
of blocklength N and rate R. More explicitly, one picks 2^ codewords as uniformly random points in 
the hypercube {0, 1}^. This means that each codeword is a string of N fair coins x^ — (xj , ■ ■ ■ , xffi) 
ioBa = l,...,2 NR . 

Once the ensemble is defined, one can estimate its average block error probability and show that it 
vanishes in the blocklength for R < C(Q). Here we will limit ourselves to providing some basic 'geometric' 
intuition of why a random code from the Shannon ensemble performs well with high probability^] 

Let us consider a particular codeword, say x^°\ and try to estimate the distance (from x^) at which 



5 The reader might notice two imprecisions with this definition. First, is not necessarily an integer: one should rather 
use [2 JV ' R ] codewords, but the difference is obviously negligible. Second, in contradiction with our definition, two codewords 
may coincide if they are independent. Again, only an exponentially small fraction of codewords will coincide and they can be 
neglected for all practical purposes. 

6 Here and in the rest of the lectures, the expression with high probability means 'with probability approaching one as 
N -» oo' 
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other codewords in £ can be found. This information is conveyed by the distance enumerator 



AC 



>(d) = # | x e €\x {0) such that d(x, x (0) ) = d } 



(2.17) 



where d(x,x') is the Hamming distance between x and x' (i.e., the number of positions in which x 
and x[ differ). The expectation of this quantity is the number of codewords different from (that is 
(2^^ — 1)) times the probability that any given codeword has distance d from x}°\ Since each entry is 
independent and different with probability 1/2, we get 



EA/>) (d) = (2 



NR 



' 2*1 d ~ 



(2.18) 



where 5 = d/N and = denotes equality to the leading exponential orderQ 

The exponent R — 1 + fj2 (<5) is plotted in Fig. 03 For 8 sufficiently small (and R < 1) this exponent 
is negative. Its first zero, to be denoted as Sqv(R), is called the Gilbert- Varshamov distance. For any 
5 < SQy(R) the expected number of codewords of distance at most NS from x} * 1 is exponentially small 
in N. It follows that the probability to find any codeword at distance smaller than N5 is exponentially 
small in N. 

Vice-versa, for d = NS, with S > 5qv(R), EACot (d) is exponentially large in N. Indeed, AC<o)(d) 
is a binomial random variable, because each of the 2^ — 1 codewords is at distance d independently 
and with the same probability. As a consequence, Af x (o)(d) is exponentially large as well with high 
probability. 

The bottom line of this discussion is that, for any given codeword x} ^ in £, the closest other codeword 
is, with high probability, at distance N(6qv(R) ± e). A sketch of this situation is provided in Fig. |U 



codewords 




Figure 4: Pictorial description of a typical code from the random code ensemble. 



Let us assume that the codeword x}°) is transmitted through a BSC(p). Denote by y € {0, 1}^ the 
channel output. By the law of large numbers d(x, y) « Np with high probability. The receiver tries to 
reconstruct the transmitted codeword from y using word MAP decoding. Using Eq. (|2.10[) . we see that 
the 'energy' of a codeword x^ (or, in more conventional terms, its log-likclihood) is given by 

N N 

E(x^) = -^logQ(^|^) = -^{l(y J =^ Q) )log(l-p)+I( ?/l ^^ a) )logp} (2.19) 
i=i i=i 

= NA(p) + 2B(p)d(x^\y), (2.20) 
7 Explicitly, we write /jv = gjv if log/jv/gjv — » 0. 
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where A(p) = —\ogp and B(p) = | log(l — p)/p. For p < 1/2, B(p) > and word MAP decoding 
amounts to finding the codeword xj-^ which is closest in Hamming distance to the channel output y. By 
the triangle inequality, the distance between y and any of the 'incorrect' codewords is > N(5gv{R) —p). 
For p < 5gv(R)/2 this is with high probability larger than the distance from x^. 

The above argument implies that, for p < S GV (R)/2, the expected block error rate of a random code 
from Shannon's ensemble vanishes as N — > oo. Notice that the channel coding theorem promises instead 
vanishing error probability whenever R < 1 — f)2(p), that is (for p < 1/2) p < Sgv(R)- The factor 2 of 
discrepancy can be recovered through a more careful argument. 

Without entering into details, it is interesting to understand the basic reason for the discrepancy 
between the Shannon Theorem and the above argument. This is related to the geometry of high dimen- 
sional spaces. Let us assume for simplicity that the minimum distance between any two codewords in € 
is at least N(5gv{R) — £)- In a given random code, this is the case for most codeword pairs. We can then 
eliminate the pairs that do not satisfy this constraint, thus modifying the code rate in a negligible way 
(this procedure is called expurgation). The resulting code will have minimum distance (the minimum 
distance among any two codewords in <£) d(<£) w N5g\j(R). 

Imagine that we use such a code to communicate through a BSC and that exactly n bits are flipped. 
By the triangular inequality, as long as n < d(€)/2, the word MAP decoder will recover the transmitted 
message for all error patterns. If on the other hand n > d(<t)/2, there are error patterns involving n bits 
such that the word-MAP decoder does not return the transmitted codeword. If for instance there exists 
a single codeword x^ at distance = 2n — 1 from x^°\ any pattern involving n out of the 2n — 1 
such that xf 1 ^ ^ will induce a decoding error. However, it might well be that most error patterns 
with the same number of errors can be corrected. 

Shannon's Theorem points out that this is indeed the case until the number of bits flipped by the 
channel is roughly equal to the minimum distance d(<t). 

3 Sparse Graph Codes 

Shannon's Theorem provides a randomized construction to find a code with 'essentially optimal' rate vs 
error probability tradeoff. In practice, however, one cannot use random codes for communications. Just 
storing the code £ requires a memory which grows exponentially in the blocklcngth. In the same vein 
the optimal decoding procedure requires an exponentially increasing effort. On the other hand, we can 
not use very short codes since their performance is not very good. To see this assume that we transmit 
over the BSC with parameter p. If the blocklcngth is N then the standard deviation of the number of 
errors contained in a block is y/ Np(l — p). Unless this quantity is very small compared to Np we have 
to either over-provision the error correcting capability of the code so as to deal with the occasionally 
large number of errors, waisting transmission rate most of the time, or we dimension the code for the 
typical case, but then we will not be able to decode when the number of errors is larger than the average. 
This means that short codes are either inefficient or unreliable (or both). 

The general strategy for tackling this problem is to introduce more structure in the code definition, 
and to hope that such structure can be exploited for the encoding and the decoding. In the next section 
we shall describe a way of introducing structure that, while preserving Shannon's idea of random codes, 
opens the way to efficient encoding/decoding. 

There are two main ingredients that make modern coding work and the two are tightly connected. 
The first important ingredient is to use codes which can be described by local constraints only. The 
second ingredient is to use a local algorithm instead of an high complexity global one (namely symbol 
MAP or word MAP decoding). In this section we describe the first component. 

3.1 Linear Codes 

One of the simplest forms of structure consists in requiring £ to be a linear subspace of {0, 1} N . One 
speaks then of a linear code. For specifying such a code it is not necessary to list all the codewords. In 
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fact, any linear space can be seen as the kernel of a matrix: 



£= {x £ {0,1} N : Hz = 0} , (3.1) 

where the matrix vector multiplication is assumed to be performed modulo 2. The matrix H is called 
the parity-check matrix. It has N columns and we let M < N denote its number of rows. Without 
loss of generality we can assume EI to have maximum rank M. As a consequence, £ is a linear space of 
dimension N — M. The rate of £ is 

(3-2) 

The a-th line in Hs — has the form (here and below © denotes modulo 2 addition) 

x ll ( a ) © • • -®x ik(a) = 0. (3.3) 

It is called a parity check. 




The parity-check matrix is conveniently represented through a factor graph (also called Tanner 
graph). This is a bipartite graph including two types of nodes: M function nodes (corresponding to the 
rows of H, or the parity-check equations) and N variable nodes (for the columns of H, or the variables). 
Edges are drawn whenever the corresponding entry in H is non- vanishing. 

Example 4: In Fig. [5] we draw the factor graph corresponding to the parity-check matrix (here 
N = 7, M = 3) 



10 10 10 1 
110 11 
1 1 1 1 



(3.4) 



In the following we shall use indices i, j, . . . for variable nodes and a,b,... for check nodes. We shall 
further denote by di (respectively, da) the set of nodes that are adjacent to variable node i (to factor 
node a). 

Remarkably, introducing the linear space structure does not deteriorate the performances of the 
resulting code. Let us introduce Shannon's parity-check ensemble: it is defined by letting the parity- 
check matrix H be a uniformly random matrix with the prescribed dimensions. Explicitly, each of 
the NM entries H a i is an independent Bernoulli random variable of mean 1/2. Probabilistic arguments 
similar to the ones for the random code ensemble can be developed for the random parity-check ensemble. 
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The conclusion is that random codes from this ensemble allow to communicate with arbitrarily small 
block error probability at any rate R < C(Q), where C(Q) is the capacity of the given BMS channel. 

Unfortunately, linearity is not sufficient to guarantee that a code admits a low-complexity decoding 
algorithm. In particular, the algorithm which we discuss in the sequel works well only for codes that can 
be represented by a sparse parity-check matrix H (i.e. a parity check matrix with O(N) non-vanishing 
entries). Notice that a given code £ has more than one representation of the form (|3.1[) . A priori one 
could hope that, given a uniformly random matrix H, a new matrix H' could be built such that H' is 
sparse and that its null space coincides with the one of EL This would provide a sparse representation of 
(C. Unfortunately, this is the case only for a vanishing fraction of matrices H, as shown by the exercise 
below. 

Exercise 2: Consider a linear code £, with blocklength N, and dimension N — M (as a linear 
space). Prove the following sequence of arguments. 

(i) The total number of binary N x M parity-check matrices is 2 . 
(it) Each code £ has 2^) JjtLi ( 2 * - l) distinct N x M parity-check matrices H. 
(Hi) The number of such matrices with at most aN non-zero entries is X^=o — 2 NMt]2 ^ ( N ~ M )), 

(iv) Conclude from the above that, for any given a, the fraction of parity-check matrices H that 
admit a sparse representation in terms of a matrix HP with at most aN ones, is of order e~ Nl 
for some 7 > 0. 

With an abuse of language in the following we shall sometimes use the term 'code' to denote a pair 
code/parity-check matrix. 

3.2 Low-Density Parity-Check Codes 

Further structure can be introduced by restricting the ensemble of parity-check matrices. Low-density 
parity-check (LDPC) codes are codes that have at least one sparse parity-check matrix. 

Rather than considering the most general case let us limit ourselves to a particularly simple family 
of LDPC ensembles, originally introduced by Robert Gallager [19 . We call them 'regular' ensembles. 
An element in this family is characterized by the blocklength N and two integer numbers k and I, with 
k > I. We shall therefore refer to it as the (k,l) regular ensemble). In order to construct a random 
Tanner graph from this ensemble, one proceeds as follows: 

1. Draw N variable nodes, each attached to I half-edges and M — Nl/k (we neglect here the possibility 
of Nl/k not being an integer) check nodes, each with k half edges. 

2. Use an arbitrary convention to label the half edges form 1 to Nl, both on the variable node side 
as well as the check node side (note that this requires that Mk = Nl). 

3. Choose a permutation 7r uniformly at random among all permutations over Nl objects, and connect 
half edges accordingly. 

Notice that the above procedure may give rise to multiple edges. Typically there will be 0(1) multiple 
edges in a graph constructed as described. These can be eliminated easily without effecting the perfor- 
mance substantially. From the analytical point of view, a simple choice consists in eliminating all the 
edges (i, a) if (i,a) occurs an even number of times, and replacing them by a single occurrence (i,a) if 
it occurs an odd number of times. 

Neglecting multiple occurrences (and the way to resolve them), the parity-check matrix corresponding 
to the graph constructed in this way does include k ones per row and I ones per column. In the sequel 
we will keep / and k fixed and consider the behavior of the ensemble as N — > 00. This implies that the 
matrix has only O(N) non-vanishing entries. The matrix is sparse. 

For practical purposes it is important to maximize the rate at which such codes enable one to 
communicate with vanishing error probability. To achieve this goal, several more complex ensembles 
have been introduced. As an example, one simple idea is to consider a generic row/column weight 
distribution (the weight being the number of non-zero elements), cf. Fig. [5] for an illustration. Such 
ensembles are usually referred to as 'irregular', and were introduced in |27j . 
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Figure 6: Factor graph of an irregular LDPC code. Variable nodes and function nodes can have any degree 
between 2 and d max . Half edges on the two sides are joined through a uniformly random permutation. 



3.3 Weight Enumerator 

As we saw in Section 12.51 the reason of the good performance of Shannon ensemble (having vanishing 
block error probability at rates arbitrarily close to the capacity), can be traced back to its minimum 
distance properties. This is indeed only a partial explanation (as we saw errors could be corrected 
well beyond half its minimum distance). It is nevertheless instructive and useful to understand the 
geometrical structure (and in particular the minimum distance properties) of typical codes from the 
LDPC ensembles defined above. 

Let us start by noticing that, for linear codes, the distance enumerator does not depend upon the 
reference codeword. This is a straightforward consequence of the observation that, for any x^ 6 £ 
the set xj ^ £ = {jrA -* ffi x : x € €} coincides with €. We are therefore led to consider the distance 
enumerator with respect to the all-zero codeword 0. This is also referred to as the weight enumerator, 

N{w) = #{x<^£ : w(x)=w], (3.5) 

where w(x) = d(x,0) is the number of non-zero entries in x. 

Let us compute the expected weight enumerator 77 (w) = E7V(w) . The final result is 

Here, F = Nl — Mk denotes the number of edges in the Tanner graph, qk(z) = ^[(l + z) k + (l — z) k ], and, 
given a polynomial p(z) and an integer n, coeff[p(z), z n ] denotes the coefficient of z n in the polynomial 
p(z). 

We shall now prove Eq. (|3.6|) . Let x £ {0, 1}^ be a binary word of length N and weight w. Notice 
that Ha; = if and only if the corresponding factor graph has the following property. Consider all 
variable nodes i such that Xi = 1, and color in red all edges incident on these nodes. Color in blue all 
the other edges. Then all the check nodes must have an even number of incident red edges. A little 
thought shows that J7(w) is the number of 'colored' factor graphs having this property, divided by the 
total number of factor graphs in the ensemble. 

A valid colored graph must have wl red edges. It can be constructed as follows. First choose w 
variable nodes. This can be done in ( ) ways. Assign to each node in this set I red sockets, and to 
each node outside the set / blue sockets. Then, for each of the M function nodes, color in red an even 
subset of its sockets in such a way that the total number of red sockets is E — wl. The number of ways 
of doing this i^l coeff [qk(z) M , z lw ]. Finally we join the variable node and check node sockets in such a 
way that colors are matched. There are (lw)\(F — lw)l such matchings out of the total number of F\ 
corresponding to different elements in the ensemble. 



This is a standard generating function calculation, and is explained in Appendix 1X1 



(3.6) 
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Figure 7: Logarithm of the expected weight enumerator for the (3, 6) ensemble in the large blocklength 
limit. Inset: small weight region. Notice that (j>(u) < for uj < lo* ~ 0.02: besides the 'all-zero' word there 
is no codeword of weight smaller than Nu>* in the code with high probability. 



Let us compute the exponential growth rate 4>(u>) of A/*(ui). Lhis is defined by 

77(w = Nlu) = e N ^ . (3.7) 

In order to estimate the leading exponential behavior of Eq. (|3.6p . we set w = Nuj and estimate the 
coeff [. ..,...] term using the Cauchy Theorem, 



^ logc7fe(z) - Iuj logz 



dz 

ss- (3 - 8) 



Here the integral runs over any path encircling the origin in the complex z plane. Evaluating the integral 
using the saddle point method we finally get J7(w) = e N< ^, where 

<j>(u) = {l-l)^{u) + -\ogq k {z)-u)l\ogz, (3.9) 
and z is a solution of the saddle point equation 

(3-10) 

k q k {z) 

The typical result of such a computation is shown in Fig. [7] As can be seen, there exists > such 
that 4>(w) < for u> E (0, a;*). This implies that a typical code from this ensemble will not have any 
codeword of weight between and iV(w* — e). By linearity the minimum distance of the code is at 
least i=a Nuj*. This implies in particular that such codes can correct any error pattern over the binary 
symmetric channel of weight w*/2 or less. 

Notice that <6(u>) is an 'annealed average', in the terminology of disordered systems. As such, it 
can be dominated by rare instances in the ensemble. On the other hand, since logAfN(Nw) — Q(N) 
is an 'extensive' quantity, we expect it to be self averaging in the language of statistical physics. In 
mathematics terms one says that it should concentrate in probability. Formally, this means that there 
exists a function $jv(w) that is non-random (i.e., does not depend upon the code) and such that 

lim F{\logj\f N {NLu)-<S> N (uj)\>N6} = 0. (3.11) 
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Further we expect that $at(w) = Ncj) q (uj) + o(N) as N — > oo. Despite being rather fundamental, both 
these statements are open conjectures. 

The coefficient </> q (w) is the growth rate of the weight enumerator for typical codes in the ensembles. 
In statistical mechanics terms, it is a 'quenched' free energy (or rather, entropy). By Jensen inequality 
^q(w) < 4>{oj). A statistical mechanics calculation reveals that the inequality is strict for general (irreg- 
ular) ensembles. On the other hand, for regular ensembles as the ones considered here, (f> q (uj) = <p(u>)- 
the annealed calculation yields the correct exponential rate. This claim has been supported rigorously 
by the results of 03 02 EI]. 

Let us finally comment on the relation between distance enumerator and the Franz-Parisi potential 
[15] . introduced in the study of glassy systems. In this context the potential is used to probe the structure 
of the Boltzmann measure. One considers a system with energy function E(x), a reference configuration 
xq and some notion of distance between configurations d{x,x'). The constrained partition function is 
then defined as 

Z(x , w) = J e~ E{x) 6{d(x , x) - w) Ax. (3.12) 

One then defines the potential ^^(uj) as the typical value of \ogZ(xo, w) when xq is a random config- 
uration with the same Boltzmann distribution and w ~ Nto. Self averaging is expected to hold here 
too: 

lim P X0 { \logZ(x ,Nuj) - $ N (w)\ > NS} = 0. (3.13) 

N—>oc 

Here N may denote the number of particles or the volume of the system and ¥ Xo { • ■ • } indicates prob- 
ability with respect to xq distributed with the Boltzmann measure for the energy function E(xq). 

It is clear that the two ideas are strictly related and can be generalized to any joint distribution of N 
variables (sci, . . . , Xjsr)- In both cases the structure of such a distribution is probed by picking a reference 
configuration and restricting the measure to its neighborhood. 

To be more specific, the weight enumerator can be seen as a special case of the Franz-Parisi potential. 
It is sufficient to take as Boltzmann distribution the uniform measure over codewords of a linear code 
(£. In other words, let the configurations be binary strings of length N, and set E(x) — if x £ £, and 
= oo otherwise. Then the restricted partition function is just the distance enumerator with respect to 
the reference codeword, which indeed does not depend on it. 



4 The Decoding Problem for Sparse Graph Codes 

As we have already seen, MAP decoding requires computing either marginals or the mode of the condi- 
tional distribution of ir_ being the channel input given output y. In the case of LDPC codes the posterior 
probability distribution factorizes according to underlying factor graph G: 

^ N M 

MCyfe) = Z (£ l) Y[Q(Vi\ x i) YlM^iia) © ■ • • © Xi k (a) = 0) • (4.1) 
\ ' y> i=l a=l 

Here (ii(a), . . . ,ik(a)) denotes the set of variable indices involved in the a-th parity check (i.e., the 
non-zero entries in the a-th row of the parity-check matrix H) . In the language of spin models, the terms 
Q{yi\xi) correspond to an external random field. The factors Jix^fa) © • • • © x i k (a) — 0) can instead be 
regarded as hard core fc-spins interactions. Under the mapping crj = ( — l) Xi , such interactions depend 
on the spins through the product <Ji 1 ( a ) • ■ ■ a i k (a)- The model (|4.ip maps therefore onto a fc-spin model 
with random field. 

For MAP decoding, minimum distance properties of the code play a crucial role in determining the 
performances. We investigated such properties in the previous section. Unfortunately, there is no known 
way of implementing MAP decoding efficiently. In this section we discuss two decoding algorithms that 
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Figure 8: Numerical simulations of bit-flipping decoding of random codes from the (5, 10) regular ensemble. 
On the left: block error rate achieved by this scheme. On the right: fraction of unsatisfied parity checks 
in the word found by the algorithm. 



exploit the sparseness of the factor graph to achieve efficient decoding. Although such strategies are sub- 
optimal with respect to word (or symbol) MAP decoding, the graphical structure can itself be optimized, 
leading to state-of-the-art performances. 

After briefly discussing bit-flipping decoding, most of this section will be devoted to message passing 
that is the approach most used in practice. Remarkably, both bit flipping as well as message passing are 
closely related to statistical mechanics. 

4.1 Bit Flipping 

For the sake of simplicity, let us assume that communication takes place over a binary symmetric channel. 
We receive the message y € {0, 1}^ and try to find the transmitted codeword x as follows: 

Bit-flipping decoder 

0. Set x(0) = y. 

1. Find a bit belonging to more unsatisfied than satisfied parity checks. 

2. If such a bit exists, flip it: Xi(t+l) = Xi(t)(Bl ■ Keep the other bits: Xj(t+1) = 
Xj (t) for all j 7^ %. 

If there is no such bit, return x(t) and halt. 

3 . Repeat steps 1 and 2 . 

The bit to be flipped is usually chosen uniformly at random among the ones satisfying the condition at 
step 1. However this is irrelevant for the analysis below. 

In order to monitor the bit-flipping algorithm, it is useful to introduce the function: 

U(t) = #{ parity-check equations not satisfied by x(t) } . (4-2) 

This is a non-negative integer, and if U(t) = the algorithm is halted and it outputs x(t). Furthermore, 
U (t) cannot be larger than the number of parity checks M and decreases (by at least one) at each cycle. 
Therefore, the algorithm complexity is O(N) (this is a commonly regarded as the ultimate goal for many 
communication problems). 

It remains to be seen if the output of the bit-flipping algorithm is related to the transmitted codeword. 
In Fig. [8] we present the results of a numerical experiment. We considered the (5, 10) regular ensemble 
and generated about 1000 random code and channel realizations for each value of the noise level p in some 
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Figure 9: Sketch of the cost function U(x) (number of unsatisfied parity checks) for a typical random 
LDPC code. Filled circles correspond to codewords, and arrows to received messages in various possible 
regimes. 



mesh. Then, we applied the above algorithm and traced the fraction of successfully decoded blocks, as 
well as the residual energy [/* = £/(i*), where t* is the total number of iterations of the algorithm. The 
data suggests that bit-flipping is able to overcome a finite noise level: it recovers the original message 
with high probability when less than about 2.5% of the bits are corrupted by the channel. Furthermore, 
the curves for the block error probability Pg f under bit-flipping decoding become steeper and steeper 
as the system size is increased. It is natural to conjecture that asymptotically, a phase transition takes 
place at a well defined noise level pbt- Pg — > for p < pbf and Pg f — > 1 for p > pbf- Numerically 
Pbf = 0.025 ± 0.005. 

This threshold can be compared with the one for word MAP decoding, that we will call p c : The 
bounds in [60] state that 0.108188 < p c < 0.109161 for the (5, 10) ensemble, while a statistical mechanics 
calculation yields p c « 0.1091. Bit-flipping is significantly sub-optimal, but it is still surprisingly good, 
given the extreme simplicity of the algorithm. 

These numerical findings can be confirmed rigorously [55 . 

Theorem 4.1. Consider a regular (I, k) LDPC ensemble and let £ be chosen uniformly at random from 
the ensemble. If I > 5 then there exists e > such that, with high probability, Bit-flipping is able to 
correct any pattern of at most Ns errors produced by a binary symmetric channel. 

Given a generic word x (i.e., a length N binary string that is not necessarily a codeword), let 
us denote, with a slight abuse of notation, by U(x) the number of parity-check equations that are not 
satisfied by x. The above result, together with the weight enumerator calculation in the previous section, 
suggests the following picture of the function U(x). If x} ^ € £, than U(x^°' ) ) = 0. Moving away from 
xj°\ U(x) will become strictly positive. However as long as d(xj°\x) is small enough, U(x) does not 
have any local minimum distinct from xj°K A greedy procedure with a starting point within such a 
Hamming radius is able to reconstruct xj-°\ As we move further away, U(x) stays positive (no other 
codewords are encountered) but local minima start to appear. Bit flipping gets trapped in such minima. 
Finally, for d(x^°\x) > Nm* new codewords, i.e., minima with U(x) = 0, are encountered. 

4.2 Message Passing 

Message-passing algorithms are iterative and have low complexity. Unlike the bit-flipping procedure in 
the previous section, the basic variables are now associated to directed edges in the factor graph. More 
precisely, for each edge (i, a) (corresponding to a non-zero entry in the parity-check matrix at row a and 
column i), we introduce two messages !/(_,„ and v a ^i- Messages are elements of some set (the message 
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Figure 10: Graphical representation of message passing updates. 



alphabet) that we shall denote by M. Depending on the specific algorithm, M can have finite cardinality, 
or be infinite, for instance M = M. Since the algorithm is iterative, it is convenient to introduce a time 
index t = 0, 1, 2, . . . and label the messages with the time at which they are updated: vf\ a and v^l^ 
(but we will sometimes drop the label below). 

The defining property of message-passing algorithms is that the message flowing from node u to v 
at a given time is a function of messages entering u from nodes w distinct from v at the previous time 
step. Formally, the algorithm is defined in terms of two sets of functions $^ a ( ■ ), *a->i( • ), that define 
the update operations at variable and function nodes as follows 

i££> = ^({^b G di\a}; yi ) , v { a l t = G 0o\t}) . (4.3) 

Notice that messages are updated in parallel and that the time counter is incremented only at variable 
nodes. Alternative scheduling schemes can be considered but we will stick to this for the sake of 
simplicity. After a pre-established number of iterations, the transmitted bits are estimated using all 
the messages incoming at the corresponding nodes. More precisely, the estimate at function i is defined 
through a new function 

xf ) (y)=M{^i i ;bedi};y i ). (4.4) 

A graphical representation of message passing updates is provided in Fig. 1101 

A specific message-passing algorithm requires the following features to be specified: 

1. The message alphabet M. 

2. The initialization {v^J, {i^J. 

3. The update functions {$i_> ( • )}, {* a ^i( ■ )}■ 

4. The final estimate functions {$;( • )}. 

The most prominent instance of a message-passing algorithm is the Belief Propagation (BP) algo- 
rithm. In this case the messages vfX a {xi) and are distributions over the bit variables Xi G {0, 1}. 
The message v£_+i(xi) is usually interpreted as the a posteriori distributions of the bit Xi given the in- 
formation coming from edge a — > i. Analogously, Vi—, a (xi) is interpreted as the a posteriori distribution 
of Xi, given all the information collected through edges distinct from (a, i). Since the messages normal- 
ization (explicitly z/j_> a (0) + i/j_> a (l) = 1) can be enforced at any time, we shall neglect overall factors 
in writing down the relation between to messages (and correspondingly, we shall use the symbol oc). 

BP messages are updated according to the following rule, whose justification we will discuss in the 
next section 

v££\xi) ex Q( yi \ Xi ) J] P^(xi), (4.5) 

b£di\a 

"aU^) <* ^ l I(xi®x jl ®~'®x ih _ l =0) J] ( 4 - 6 ) 

{xj} jeda\i 
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Figure 11: Factor graph of a regular LDPC code, and notation for the belief propagation messages. 



where we used ■ ■ ■ ,jk—i) to denote the neighborhood da of factor node a. After any number of 

iterations the single bit marginals can be estimated as follows 



vf +1 \ Xi ) <x OWlIto. (4.7) 

bedi 

The corresponding MAP decision for bit i (sometimes called 'hard decision', while Vi(xi) is the 'soft 
decision') is 

Xt — argmax vf\xi) . (4.8) 

Xi 

Notice that the above prescription is ill-defined when ^(0) = It turns out that it is not really 

important which rule to use in this case. To preserve the 0—1 symmetry, we shall assume that the 
decoder returns = or = 1 with equal probability. 

Finally, as initial condition one usually takes z?^_ > 1 i (-) to be the uniform distribution over {0,1} 
(explicitly i£3 (0) = ££3 (1) = 1/2). _ 

Since for BP the messages arc distributions over binary valued variables, they can be described by a 
single real number, that is often chosen to be the bit log-likelihoodH 

^ a= 2 l0g ^-(I)' M -=2 l0g p-(I)- ^ 

We refer to Fig. QT] for a pictorial representation of these notations. We further introduce the channel 
log-likelihoods 

The BP update equations (|4. 5|) . (|4.6[) read in this notation 



E «W i = atauh{ J] tanh^l a }. (4.11) 



b£di\a j£da\i 



In this language the standard message initialization would be u„_Jj = 0. Finally, the overall log-likelihood 



at bit i is obtained by combining all the incoming messages in agreement with Eq. (|4.7|) . One thus gets 
the decision rule 



I o if^ + s^A >o, 



9 The conventional definition of log- likelihoods does not include the factor 1/2. We introduce this factor here for uniformity 
with the statistical mechanics convention (the h's and it's being analogous to effective magnetic fields). 
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Notice that we did not commit to any special decision if B t + J2 bGdi u b ^ i = 0. To keep complete 
symmetry we'll establish that the decoder returns or 1 with equal probability in this case. 

4.3 Correctness of Belief Propagation on Trees 

The justification for the BP update equations (|4.5[) . (|4.6[) lies in the observation that, whenever the 
underlying factor graph is a tree, the estimated marginal v\ (x^) converges after a finite number of 
iterations to the correct one fii(xi). In particular, under the tree assumption, and for any t sufficiently 
large, (y) coincides with the symbol MAP decision. 

In order to prove this statement, consider a tree factor graph G. Given a couple of adjacent nodes 
u,v, denote by G(u — > v) the subtree rooted at the directed edge u — > v (this contains all that can be 
reached from v through a non-reversing path whose first step is v — > u). If i is a variable index and a 
a parity-check index, let Hi-* a { ■ ) be the measure over x = {xj : j € G(i — > a)}, that is obtained by 
retaining in Eq. (|4.1[) only those terms that are related to nodes in G(i — > a): 

fii^a{x) = JJ Q(yi\Xi) Yl l (xii(b) ® ■ ■ ■ ® x ik{b) =0). (4.13) 

^ ^ a ' jeG(i^a) b£G(i->a) 

The measure /Lt a -»i( ' ) is defined analogously for the subtree G{a — > i). The marginals Hi^ a {xi) (respec- 
tively fi a ^i{xi)) are easily seen to satisfy the recursions 

Hi^ a {xi) cx Q{yi\xi) JJ_ fib^i(xi), (4-14) 

b£di\a 

) CX ^ I(Xi 8 X h ® ■ ■ ■ © X J - fc _ 1 = 0) J| fXj^aiXj), (4.15) 
{xj} jeda\i 

which coincide, apart from the time index, with the BP recursion (|4.5p . (|4.6[) . That such recursions 
converges to {fj,i^ a (xi), jLt a _>i(^)} follows by induction over the tree depth. 

In statistical mechanics equations similar to (14. 14[) . (|4.15[) are often written as recursions on the 
constrained partition function. They allow to solve exactly models on trees. However they have been 
often applied as mean-field approximation to statistical models on non-tree graphs. This is often referred 
to as the Bethe-Peierls approximation [5]. 

The Bethe approximation presents several advantages with respect to 'naive-mean field' [ST] (that 
amounts to writing 'self-consistency' equations for expectations over single degrees of freedom). It 
retains correlations among degrees of freedom that interact directly, and is exact on some non-empty 
graph (trees). It is often asymptotically (in the large size limit) exact on locally tree- like graphs. Finally, 
it is quantitatively more accurate for non-tree like graphs and offers a much richer modeling palette. 

Within the theory of disordered systems (especially, glass models on sparse random graphs), Eqs. (|4. 14[) 
and (I4.15P are also referred to as the cavity equations. With respect to Bethe-Peierls, the cavity ap- 
proach includes a hierarchy of ('replica symmetry breaking') refinements of such equations that aim at 
capturing long range correlations [37) . This will briefly described in Section [5j 

We should finally mention that several improvements over Bethe approximation have been developed 
within statistical physics. Among them, Kikuchi's cluster variational method [53] is worth mentioning 
since it motivated the development of a 'generalized belief propagation' algorithm, which spurred a lot 
of interest within the artificial intelligence community [61] . 

4.4 Density Evolution 

Although BP converges to the exact marginals on tree graphs, this says little about its performances 
on practical codes such as the LDPC ensembles introduced in Section [3] Fortunately, a rather precise 
picture on the performance of LDPC ensembles can be derived in the large blocklength limit N — > oo. 
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The basic reason for this is that the corresponding random factor graph is locally tree-like with high 
probability if we consider large blocklengths. 

Before elaborating on this point, notice that the performance under BP decoding (e.g., the bit error 
rate) is independent on the transmitted codeword. For the sake of analysis, we shall hereafter assume 
that the all-zero codeword has been transmitted. 

Consider a factor graph G and let (i, a) be one of its edges. Consider the message z^-L sen t by the 
BP decoder in iteration t along edge (i,a). A considerable amount of information is contained in the 
distribution of with respect to the channel realization, as well as in the analogous distribution for 

^a-Li- To see this, note that under the all-zero codeword assumption, the bit error rate after t iterations 
is given by 

i n 

p b ) = -E p {^^-i; 6G9 *^)^ } ■ ( 4 - 16 ) 

i=l 

Therefore, if the messages vj^A are independent, then the bit error probability is determined by the 
distribution of z?„_,j. 

Rather than considering one particular graph (code) and a specific edge, it is much simpler to take 
the average over all edges and all graph realizations. We thus consider the distribution a[ N \ ■ ) of vf\ a 
with respect to the channel, the edges, and the graph realization. While this is still a quite difficult 
object to study rigorously, it is on the other hand possible to characterize its large blocklength limit 
a f ( ■ ) = limjv a* ( ■ )• This distribution satisfies a simple recursion. 

It is convenient to introduce the directed neighborhood of radius r of the directed edge i — > a in G, 
call it Bj_, a (r; G). This is defined as the subgraph of F that includes all the variable nodes that can be 
reached from i through a non-reversing path of length at most r, whose first step is not the edge (i, a). 
It includes as well all the function nodes connected only to the above specified variable nodes. In Fig. [T3] 
we reproduce an example of a directed neighborhood of radius r — 3 (for illustrative purposes we also 
include the edge (i, a)) in a (2, 3) regular code. 

If F is the factor graph of a random code from the (k,l) LDPC ensemble, then Bj_> (r; F) is with 
high probability a depth-r regular tree with degree I at variable nodes and degree k at check nodes (as 
in Fig. 1131 where / = 2 and k = 3). The basic reason for this phenomenon is rather straightforward. 
Imagine to explore the neighborhood progressively, moving away from the root, in a breadth first fashion. 
At any finite radius r, about c r /N vertices have been visited (here c = (k — 1)(Z — 1)). The vertices 
encountered at the next layer will be 'more or less' uniformly random among all the ones not visited so 
far. As a consequence they will be distinct with high probability, and Bi_> a (r + 1; G) will be a tree as 
well. This argument breaks down when the probability that two of the 0(c r ) new vertices coincide, that 
is for c 2r = 9(iV)@ This is equivalent to r ~ \ log c N. 

The skeptical reader is invited to solve the following exercise. 

Exercise 3: In order to illustrate the above statement, consider the example of a random code from 
the regular (2,3) ensemble (each variable has degree 2 and each check has degree 3). The three 
possible radius- 1 neighborhoods appearing in the associated factor graph are depicted in Fig. [Ml 

(a) Show that the probability that a given edge (i, a) has neighborhoods as in (B) or (C) is 
0(1/N). 

(b) What changes for a generic radius r? 

For illustrative reasons, we shall occasionally add a 'root edge' to Bi_> a (r; G), as for i — ► a in Fig. [T3l 

Now consider the message vf\ a . This is a function of the factor graph G and of the received message 
y. However, a moment's thought shows that it will depend on G only through its directed neighborhood 
Bi^ a (t + 1; G), and only on the received symbols yj, j € B^_, (t; G). 



This is the famous birthday problem. The probability that two out of a party of n peoples were born on the same day of 
the year, scales like n 2 /N for n 2 <C N (N is here the number of days in a year). 
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Figure 12: Evolution of the probability density functions of an u^ t+1 ^ for an irregular LDPC code used 
over a gaussian channel. From top to bottom t = 0, 5, 10, 50, and 140. 
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In view of the above discussion, let us consider the case in which Bj_> a (t+1; G) is a (k, £)-regular tree. 
We further assume that the received symbols yj are i.i.d. with distribution Q(y\0), and that the update 
rules (|4.3[) do not depend on the edge we are considering (i.e., $j_> a ( ' ) = ) an d ^i-> a ( • ) = ■ ) 
independent of i, a). 

Let v {t ^ be the message passed through the root edge of such a tree after t BP iterations. Since the 

actual neighborhood Bj_ >a (t+ 1; G) is with high probability a tree, — ► as iV — * oo. The symbol 

-i denotes convergence in distribution. In other words, for large blocklengths, the message distribution 
after t iterations is asymptotically the same that we would have obtained if the graph were a tree. 

Consider now a (k, Z)-regular tree, and let j — > b an edge directed towards the root, at distance d 
from it. It is not hard to realize that the message passed through it after r — d — 1 (or more) iterations is 
distributed as i/( r ~ d ~ 1 ) . Furthermore, if j% — > 61 and j% — » hi are both directed upwards and none belongs 
to the subtree rooted at the other one, then the corresponding messages are independent. Together with 
Eq. (|4.3|) . these observation imply that 

4 . . .^ v y) , P<«> 4 . . . , ^) . (4.17) 

Here , . . . ,v[ t } 1 are i.i.d. copies of and uf 1 , . . . , v^_-y i.i.d. copies of v^ 1 . Finally, y is a received 
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symbol independent from the previous variables and distributed according to Q(y\0). 

Equations (|4.17[) , or the sequence of distributions that they define, are usually referred to as density 
evolution. The name is motivated by the identification of the random variables with their densities 
(even if these do not necessarily exist). They should be parsed as follows (we are refer here to the first 
equation in (|4.17p : an analogous phrasing holds for the second): pick I — 1 i.i.d. copies and y with 

distribution Q(y\0), compute <&(y^\ . . . , P;_\; y)- The resulting quantity will have distribution v( t+1 \ 
Because of this description, they are also called 'recursive distributional equations'. 

Until this point we considered a generic message passing procedure. If we specialize to BP decoding, 
we can use the parametrization of messages in terms of log-likelihood ratios, cf. Eq. (|4.9|1 . and use the 
above arguments to characterize the limit random variables and 14W. The update rules (|4.1ip then 
imply 

= B + + • • • + u[ t } 1 , u<*> =atanh{tanh/ii t) ---tanh4-i} ■ ( 4 - 18 ) 

Here , . . . , are i.i.d. copies of . . . , /i^-i are i-i-d. copies of and B = | log < q$\X ; 

where y is independently distributed according to Q(y\Q). It is understood that the recursion is initiated 
with = 0. 

Physicists often write distributional recursions explicitly in terms of densities. For instance, the first 
of the equations above reads 

a t +i(h)= m da t (u 6 ) dp(B) 5 (h - B - £ u b J , (4.19) 

^ 6=1 V 6=1 / 

where a*t( • ) denotes the density of , and p( • ) the density of B. We refer to Fig. rT2]for an illustration 
of how the densities a t ( • ), a t ( • ) evolve during the decoding process. 

In order to stress the importance of density evolution notice that, for any continuous function /(#), 



1 N 

jv lim ) E{-X;/(^ )} =n.f(h (t) )}, (4.20) 



where the expectation is taken with respect to the code ensemble. Similar expressions can be obtained 
for functions of several messages (and are particularly simple when such message are asymptotically 
independent). In particular^, if we let p^'*-* be the expected (over an LDPC ensemble) bit error rate 
for the decoding rule (|4.12[) . and let = lhriAr^oo p^'*- 1 be its large blocklength limit. Then 



P^ = P{B + h?> + ■■■ + hf } < 0} + - ¥{B + h { *> + ■■■ + hf> = 0} , (4.21) 
where h\ , . . . ,h\ ' are i.i.d. copies of fcW. 
4.5 The Belief Propagation Threshold 

Density evolution would not be such an useful tool if it could not be simulated efficiently. The idea 
is to estimate numerically the distributions of the density evolution variables {h^\u^}. As already 
discussed this gives access to a number of statistics on BP decoding, such as the bit error rate Pj^ after 
t iterations in the large blocklength limit. 

A possible approach consists in representing the distributions by samples of some fixed size. Within 
statistical physics this is sometimes called the population dynamics algorithm (and made its first ap- 
pearance in the study of the localization transition on Cay ley trees [46]). Although there exist more 



11 The suspicious reader will notice that this is not exactly a particular case of the previous statement, because f(x) = I(x < 
0) + 2^{ x — 0) i s n °t a continuous function. 
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Figure 15: The performances of two LDPC ensembles as predicted by a numerical implementation of 
density evolution. On the left, the (3, 6) regular ensemble. On the right, an optimized irregular ensemble. 
Dotted curves refer (from top to bottom) to t = 0, 1, 2, 5, 10, 20, 50 iterations, and bold continuous lines 
to the limit t — > oo. In the inset we plot the expected conditional entropy KH(Xi\uf^). 



efficient alternatives in the coding context (mainly based on Fourier transform, see [49J, [48] ) , we shall 
describe population dynamics because it is easily programmed. 

Let us describe the algorithm within the setting of a general message passing decoder, cf. Eq. (|4.17|) . 
Given an integer Af 3> 1, one represent the messages distributions with two samples of size Af: = 
{vx , ■ ■ ■ , and 5)3 W = {V^' , . . . ,Pw}. Such samples are used as proxy for the corresponding 

distributions. For instance, one would approximate an expectation as 

AT 

E/(^))«^^/(^). (4.22) 

i— 1 

The populations are updated iteratively. For instance <p( t+1 ) is obtained from 5)}(*) by generating 
Vi 1 ' , . . . , independently as follows. For each i 6 [A/], draw indices b\{i), . . . ,bi(i) indepen- 

dently and uniformly at random from [TV], and generate iji with distribution Q(y\0). Then compute 
vf +r > = $({P 6 ( * ) (i) };j/i) and store it in ^ t+1 l 

An equivalent description consists in saying that we proceed as if exactly represents the distri- 
bution of 14W (which in this case would be discrete). If this was the case, the distribution of h^ t+1 ^ would 
be composed of |^4| • TV' -1 Dirac deltas. In order not to overflow memory, the algorithm samples Af 
values from such a distribution. Empirically, estimates of the form (|4.22p obtained through population 
dynamics have systematic errors of order Af" 1 and statistical errors of order A/" -1 / 2 with respect to the 
exact value. 

In Fig.[T5]we report the results of population dynamics simulations for two different LDPC ensembles, 
with respect to the BSC. We consider two performance measures: the bit error rate and the bit 
conditional entropy H^K The latter is defined as 

1 N 

H (t) = £m ^EffpQlFf), (4.23) 
1=1 

and encodes the uncertainty about bit Xi after t BP iterations. It is intuitively clear that, as the algorithm 
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k 


R 


Pd 


Shannon limit 


3 


4 


1/4 


0.1669(2) 


0.2145018 


3 


5 


2/5 


0.1138(2) 


0.1461024 


3 


6 


1/2 


0.0840(2) 


0.1100279 


4 


6 


1/3 


0.1169(2) 


0.1739524 



Table 1: Belief propagation thresholds for a few regular LDPC ensembles. 



progresses, the bit estimates improve and therefore and should be monotonically decreasing 
functions of the number of iterations. Further, they are expected to be monotonically increasing functions 
of the crossover probability p. Both statement can be easily checked on the above plots, and can be 
proved rigorously as well. 

Since is non-negative and decreasing in t, it has a finite limit 

Pg p = lim , (4.24) 

t — >oo 

which is itself non-decreasing in p. The limit curve P b p is estimated in Fig. [T5] by choosing t large 

enough so that P b '^ is independent of t within the numerical accuracy. 

Since P pp = P pp (p) is a non-decreasing function of p, one can define the BP threshold 

p d = sup { p e [0, 1/2] : P pp (p) = } . (4.25) 

Analogous definitions can be provided for other channel families such as the BEC(e). In general, the 
definition (|4.25[) can be extended to any family of BMS channels BMS(p) indexed by a real parameter 
p 6 J, J C R being an interval (obviously the sup will be then taken over p G I). The only condition 
is that the family is 'ordered by physical degradation'. We shall not describe this concept formally, but 
limit ourselves to say that that p should be an 'honest' noise parameter, in the sense that the channel 
worsen as p increases. 

Analytical upper and lower bounds can be derived for pd- In particular it can be shown that it is 
strictly larger than (and smaller than 1/2) for all LDPC ensembles with minimum variable degree at 
least 2. Numerical simulation of density evolution allows to determine it numerically with good accuracy. 
In Table [4751 we report the results of a few such results. 

Let us stress that the threshold pd has an important practical meaning. For any p < Pd one can 
achieve arbitrarily small bit error rate with high probability by just picking one random code from the 
ensemble LDPC and using BP decoding and running it for a large enough (but independent of the 
blocklength) number of iterations. For p > Pd the bit error rate is asymptotically lower bounded by 
P pp (p) > for any fixed number of iterations. In principle it could be that after, let's say n a , a > 
iterations a lower bit error rate is achieved. However simulations show quite convincingly that this is 
not the case. 

In physics terms the algorithm undergoes a phase transition at pd- At first sight, such a phase 
transition may look entirely dependent on the algorithm definition and not 'universal' in any sense. As 
we will discuss in the next section, this is not the case. The phase transition at pd is somehow intrinsic to 
the underlying measure n(x), and has a well studied counterpart in the theory of mean field disordered 
spin models. 

Apart from the particular channel family, the BP threshold depends on the particular code ensemble, 
i.e. (for the case considered here) on the code ensemble. It constitutes therefore a primary measure 
of the 'goodness' of such a pair. Given a certain design rate R, one would like to make pd as large as 
possible. This has motivated the introduction of code ensembles that generalize the regular ones studied 
here (starting from 'irregular' ones). Optimized ensembles have been shown to allow for exceptionally 
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good performances. In the case of the erasure channel, they allowed to saturate Shannon's fundamental 
limit [27] • This is an important approach to the design of LDPC ensembles. 

Let us finally mention that the BP threshold was defined in Eq. (I4.25|) in terms of the bit error rate. 
One may wonder whether a different performance parameter may yield a different threshold. As long as 
such parameter can be written in the form -L f(h\^) this is not the case. More precisely 

p d = sup {pel : h {t) -i +00 I , (4.26) 

where, for the sake of generality we assumed the noise parameter to belong to an interval ICR. In 
other words, for any p < p^ the distribution of BP messages becomes a delta at plus infinity. 

4.6 Belief Propagation versus MAP Decoding 

So far we have seen that detailed predictions can be obtained for the performance of LDPC ensembles 
under message passing decoding (at least in the large blocklength limit). In particular the threshold 
noise for reliable communication is determined in terms of a distributional recursion (density evolution) . 
This recursion can in turn be efficiently approximated numerically, leading to accurate predictions for 
the threshold. 

It would be interesting to compare such predictions with the performances under optimal decoding 
strategies. Throughout this section we shall focus on symbol MAP decoding, which minimizes the bit 
error rate, and consider a generic channel family {BMS(p)} ordered^ by the noise parameter p. 

Given an LDPC ensemble, let Pj, be * ne expected bit error rate when the blocklength is N. The 
MAP threshold p c for such an ensemble can be defined as the largest (or, more precisely, the supremum) 
value of p such that lirajv— Pl^ = 0. In other words, for any p < p c one can communicate with 
an arbitrarily small error probability, by using a random code from the ensemble, provided N is large 
enough. 

By the optimality of MAP decoding, pd < p c . In coding theory some techniques have been developed 
to prove upper and lower bounds on p c [TIB [52] . In particular it is easy to find ensembles for which there 
exist a gap between the two thresholds (namely pd < p c strictly). Consider for instance (k,l) regular 
ensembles with a fixed ratio l/k = 1 — R. It is then possible to show that, as k, I — > 00, the BP threshold 
goes to while the MAP threshold approaches the Shannon limit. 

This situation is somewhat unsatisfactory. The techniques used to estimate pd and p c are completely 
different. This is puzzling since the two thresholds can be extremely close and even coincide for some 
ensembles. Furthermore, we know that pd < Pc by a general argument (optimality of MAP decoding), 
but this inequality is not 'built in' the corresponding derivations. Finally, it would be interesting to have 
a sharp estimate for p c . 

It turns out that a sharp characterization of p c can be obtained through statistical mechanics tech- 
niques 42, 58, 38 . The statistical mechanics result has been proved to be a rigorous upper bound for 
general code ensembles, and it is conjectured to be tight pM 135] . 

The starting point is to consider the conditional entropy of the channel input x given the output y, 
Hn(X\Y). As shown in Eq. (|2.f 6[) this is given by the expectation of the log partition function appearing 
in Eq. (|4.1|) (apart from a trivial additive factor). 

Let fjv = EHn(X\Y)/N denote the entropy density averaged over the code ensemble. Intuitively 
speaking, this quantity allows to estimate the typical number of inputs with non-negligible probability 
for a given channel output. If f/v is bounded away from as N — > 00, the typical channel output 
corresponds to an exponential number of (approximately) equally likely inputs. If on the other hand 
f/v — > 0, the correct input has to be searched among a sub-exponential number of candidates. This leads 
us to identif\P^I the MAP threshold as the largest noise level such that fjv — > as N — > 00. 



Such that the channel worsen as p increases. Examples are the binary symmetric or binary erasure channels. 
A rigorous justification of this identification can be obtained using Fano's inequality. 
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The Bethe free energy provides a natural way to approximate log-partition functions on sparse graphs. 
It is known to be exact if the underlying graph is a tree and its stationary points are in correspondence 
with the fixed points of BP. In statistical physics terms, it is the correct variational formulation for the 
Bethe Peierls approximation. In random systems which are locally tree like, it is normally thought to 
provide the correct N — > oo limit unless long range correlations set in. These are in turn described 
through 'replica symmetry breaking' (see below). 

As many mean field approximations, the Bethe approximation can be thought of as a way of writing 
the free energy as a function of a few correlation functions. More specifically, one considers the single- 
variable marginals {bi(xi) : i S {1, ...,iV}}, and the joint distributions of variables involved in a 
common check node {b a (x a ) : a g {1, . . . , M}}. In the present case the Bethe free energy reads 



JV 



F B (b) = -^53& i (:E i )logQ(y i |:E i ) + 

i— 1 Xi 

M N 

+ J2^2b a (x a )\ogb a (x a ) - J2(\di\ - l)J2b t (xi)\og 2 h(xi) 



(4.27) 



The marginals {&;( • )}, {b a ( ■ )} are regarded as variables. They are constrained to be probability distri- 
butions (hence non-negative) and to satisfy the marginalization conditions 



E 

Xj, j£da\i 



b a {x a ) = bi(xi) Vi € da. 



}^bj(xj) = 1 



Vi. 



(4.28) 



Further, in order to fulfill the parity-check constraints b a (x a ) must be forced to vanish unless x^m © 
' ' ' © x ia(k) = (as usual we use the convention OlogO = 0). Since they do not necessarily coincide with 
the actual marginals of /i( • ), the {b a }, {b{\ are sometimes called beliefs. 

Approximating the log-partition function — log Z(y) requires minimizing the Bethe free energy FbQ)) . 
The constraints can be resolved by introducing Lagrange multipliers, that are in turn expressed in terms 
of two families of real valued messages u = {u a _>j}, h = {hi^ a }. If we denote by P u (x) the distribution 
of a bit x whose log likelihood ratio is u (in other words P u (0) 
the resulting beliefs read 



l/(l + e^"), P tt (l) = e-^/(l + e-^)) 



b a (x a ) = — lafe) IT Ph^aixj) i bi(xi) = —Q(yi\xi) Y\ Pu a ^A x i) 



(4.29) 



jeda 



a G Oi 



where we introduced the shorthand I Q (x) to denote the indicator function for the a-th parity check being 
satisfied. Using the marginalization conditions (|4.28[) as well as the stationarity of the Bethe free energy 
with respect to variations in the beliefs, one obtains the fixed point BP equations 



hi — *n Pi 



b£di\a 



atanh < tanh hj^ a > . 

j£da\i 



(4.30) 



These in turn coincide for with the fixed point conditions for belief propagation, cf. Eqs. (|4.1ip . 

The Bethe free energy can be written as a function of the messages by plugging the expressions (I4.29P 
into Eq. (14.271) . Using the fixed point equations, we get 



N 



8=1 



^p Ua ^M)Ph^ a {pi) 
22,Q{vi\xi) Y[ p u a -,,(xi 



(4.31) 



a=l 



y.ux) n p^ji 



27 



We are interested in the expectation of this quantity with respect to the code and channel realization, 
in the N — ► oo limit. We assume that messages are asymptotically identically distributed u a —>i = 
u, hi^ a = h, and that messages incoming in the same node along distinct edges are asymptotically 
independent. Under these hypotheses we get the limit 



hm i-E f b (u, u) = -<t> u , h + J2 Q(y\°) lo s 2 Q(y\o) 

N—^oo IV — 



(4.32) 



where 



-lE Uth log 2 

-r E {^} lo g2 



^P u (a;)P /l (a;) 



■ E y E {Ut} log 2 



i=l 



(4.33) 



Notice that the random variables u, h are constrained by Eq. (|4. 30(1 . which must be fulfilled in dis- 
tributional sense. In other words u, h must form a fixed point of the density evolution recursion (|4.18[) . 
Given this proviso, if the above assumptions are correct and the Bethe free energy is a good approxima- 
tion for the log partition function one expects the conditional entropy per bit to be limjv^oo fjv = <fi u ,h- 
This guess is supported by the following rigorous result. 

Theorem 4.2. If u, h are symmetric random variables satisfying the distributional identity u = 
atanh | HiLi ta- n h hi | , then 



lim fjy > (f) u .h 



(4.34) 

It is natural to conjecture that the correct limit is obtained by optimizing the above lower bound, 



lim fjy = sup 

N-+oo uh 



(4.35) 



where, once again the sup is taken over the couples of symmetric random variables satisfying u = 
atanh tanh/ij|. In fact it is easy to show that, on the fixed point, the distributional equation 

h = B + y^j— i u a must be satisfied as well. In other words the couple u, h must be a density evolution 
fixed point. 

This conjecture has indeed been proved in the case of communication over the binary erasure channel 
for a large class of LDPC ensembles (including, for instance, regular ones). 

The expression (|4.35|) is interesting because it bridges the analysis of BP and MAP decoding. For 
instance, it is immediate to show that it implies pd < p c - 

Exercise 4: This exercise aims at proving the last statement. 

(a) Recall that u,h = +oo constitute a density evolution fixed point for any noise level. Show 
that <fih.u — on such a fixed point. 

(b) Assume that, if any other fixed point exists, then density evolution converges to it (this can 
indeed be proved in great generality). 

(c) Deduce that pd < p c ■ 

Evaluating the expression (|4.35|) implies an a priori infinite dimensional optimization problem. In 
practice good approximations can be obtained through the following procedure: 

1. Initialize /i, u to a couple of symmetric random variables h^°\ u^°\ 
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Shannon limit 
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4 


1/4 


0.2101(1) 


0.2145018 
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5 


2/5 


0.1384(1) 


0.1461024 
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6 


1/2 


0.1010(2) 


0.1100279 


4 


6 


1/3 


0.1726(1) 


0.1739524 



Table 2: MAP thresholds for a few regular LDPC ensembles and communication over the BSC(p). 



2. Implement numerically the density evolution recursion (|4.18p and iterate it until an approximate 
fixed point is attained. 

3. Evaluate the functional (f> u ,h on such a fixed point, after enforcing u = atanh |n<=i tanh/ii| 
exactly. 

The above procedure can be repeated for several different initializations u^°\ h^°\ The largest of the 
corresponding values of 4>u,h is then picked as an estimate for limjv^oo fjv- 

While his procedure is not guaranteed to exhaust all the possible density evolution fixed points, it 
allows to compute a sequence of lower bounds to the conditional entropy density. Further, one expects 
a small finite number of density evolution fixed points. In particular, for regular ensembles and p > pd, 
a unique (stable) fixed point is expected to exist apart from the no-error one u,h= +oo. In Table [4~6l 
we present the corresponding MAP thresholds for a few regular ensembles. 

For further details on these results, and complete proofs, we refer to |39j . Here we limit ourselves to 
a brief discussion why the conjecture (|4.35p is expected to hold from a statistical physics point of view. 

The expression (|4.35[) corresponds to the 'replica symmetric ansatz' from the present problem. This 
usually breaks down if some form of long-range correlation ('replica symmetry breaking') arises in the 
measure fi(-). This phenomenon is however not expected to happen in the case at hand. The technical 
reason is that the so-called Nishimori condition holds for /i( • ) [35]. This condition generally holds for a 
large family of problems arising in communications and statistical inference. While Nishimori condition 
does not provide an easy proof of the conjecture (|4.35[) . it implies a series of structural properties of 
/i( • ) that are commonly regarded as incompatible with replica symmetry breaking. 

Replica symmetry breaking is instead necessary to describe the structure of 'metastable states' [T7] . 
This can be loosely described as very deep local minima in the energy landscape introduced in Section al! 
Here 'very deep' means that Q(N) bit flips are necessary to lower the energy (number of unsatisfied 
parity checks) when starting from such minima. As the noise level increases, such local minima become 
relevant at the so called 'dynamic phase transition'. It turns out that the critical noise for this phase 
transition coincides with the BP threshold pd- In other words the double phase transition at pd and p c 
is completely analogous to what happens in the mean field theory of structural glasses (see for instance 
Parisi's lectures at this School). Furthermore, this indicates that pd has a 'structural' rather than purely 
algorithmic meaning. 

5 Belief Propagation Beyond Coding Theory 

The success of belief propagation as an iterative decoding procedure has spurred a lot of interest in its 
application to other statistical inference tasks. 

A simple formalization for this family of problems is provided by factor graphs. One is given a factor 
graph G = (V, F, E) with variable nodes V, function nodes F, and edges E and considers probability 
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Figure 16: Factor graph representation of the satisfiability formula (|5.5p . Circle correspond to variables 
and squares to clauses. Edges are represented as dashed line if the variable is negated in the corresponding 
clause. 



distributions that factorize accordingly 

aeF 

Here the variables Xi take values in a generic finite alphabet X, and the compatibility functions ip a : 
X da — > M + encode dependencies among them. The prototypical problem consists in computing marginals 
of the distribution /i( ■ ), e.g., 

Hi(xi) = ^2fi(x). (5.2) 

Belief propagation can be used to accomplish this task in a fast and distributed (but not necessarily 
accurate) fashion. The general update rules read 

^(^oc n 3&o*). ^i(^)«E^^a) n ^i(^)' ( 5 - 3 ) 

b£di\a jGda\i 

Messages are then used to estimate local marginals as follows 

The basic theoretical question is of course to establish a relation, if any between • ) and Fj(- ). 

As an example, we shall consider satisfiability [23j . Given iV Boolean variables ajj, i G {1, . . . ,N}, 
Xi G {True, False}, a formula is the logical expression obtained by taking the AND of M clauses. Each 
clause is the logical OR of a subset of the variables or their negations. As an example, consider the 
formula (here x~i denotes the negation of Xi) 

(x\ V X2 V Xi)A(xi V X2)/\{X2 Vl4 V X5)A(xi V X2 V ^5)A(xi V X3 V X5) . (5-5) 

An assignment of the N variables satisfies the formula if, for each of the clause, at least one of the 
involved literals (i.e. either the variables or their negations) evaluates to True. 

A satisfiability formula admits a natural factor graph representation where each clause is associated 
to a factor node and each variable to a variable node. An example is shown in Fig. 1161 Admitting that 
the formula has at least one satisfying assignment, it is natural to associate to a J 7 the uniform measure 



(5.4) 
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over such assignments fi^(x). It is easy to realize that such a distribution takes the form (|5.ip where 
the compatibility function ^p a {xda) takes value 1 if the assignment x satisfies clause a and otherwise. 

Satisfiability, i.e., the problem of finding a solution to a satisfiability formula or proving that it is 
usatisfiable, is one of the prototypical NP-complete problems. Given a satisfiability formula, computing 
marginals with respect to the associated distribution /^jf( ■ ) is relevant for tasks such as counting the 
number of solutions of T or sampling them uniformly. These are well known ff-P complete problems 

The currently best algorithms for solving random instances for the i^-SAT problem are based on a 
variant of BP, which is called survey propagation I33[ lllj . 

5.1 Proving Correctness through Correlation Decay 

A simple trick to bound the error incurred by BP consists in using the correlation decay properties 
[Tl [57] of the measure /z( • ). Let i £ {1, . . . , N} be a variable index and denote by x t the set of variables 
sitting on nodes at distance t from i. Further, denote by x >t the set of variables whose distance from i 
is at least t. Then the local structure of the probability distribution (|5.ip 

Hi{xi) =^2fJ,(x i \x> t )fi(x> t ) = J2fi(xi\x t )fi(x t ). (5.6) 

Let Bi(t) denote the subgraph induced by nodes whose distance from i is at most t, and dBi(t) its 
boundary (nodes whose distance from i is exactly t). Further, for any j G dBi(t) let a(j) be the unique 
function nodes inside Bi(t) that is adjacent to j. It is intuitively clear that belief propagation computes 
the marginal at i as if the graph did not extend beyond (t) . More precisely, if the initial condition 
v^l a (xi) is properly normalized, then we have the exact expression 

As a consequence of Eq. (|5.6[) and (|5.7[) we have 

\Hi{xi) -vf'{xi)\ < sup \fi(xi\x t ) - ^(xi\x! t )\ . (5.8) 

This provides an upper bound on the error incurred by BP when computing the marginal of Xi base on 
the local structure of the underlying graph in terms of the influence of far away variables. To make things 
fully explicit, assume that the graph has girtfi^g and that sup x ^ \/i(xi\x t ) — f^(xi\x^ t )\ < exp(— nt) for 
some positive k. This implies 

\lH(xi)-V?\xi)\ <e^/\ (5.9) 

As an example of such error estimates, we shall consider random k- satisfiability [16] . This is a 
standard model to generate 'synthetic' satisfiability formulae. It amounts to picking a formula uniformly 
at random among all the ones including N variables and M = Na fc-clauses (a fc-clause is a clause that 
involve exactly k distinct variables). We shall of course limit to k > 2, the case k = 1 being trivial. 

Consider a uniformly random variable node in the factor graph associated to a random formula, and 
its depth-t neighborhood Bj(t). Proceeding as in the previous section it is not hard to show that, for 
any fixed t, Bi(t) is with high probability (as N — > oo) a tree. An appropriate model the distribution 
of such a tree, is given by the tree ensemble T*(f) described as follows. For t = 0, it is the graph 
containing a unique variable node. For any t > 1, start by a single variable node (the root) and add 

Z = Poi sson(fca) clauses, each one including the root, and k—1 new variables (first generation variables). 



The notation #-P refers to the hardness classification for counting problems. 
Recall that the girth of a graph is the length of its shortest cycle. 
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For each one of the I clauses, the corresponding literals are non-negated or negated independently with 
equal probability. If t > 2, generate an independent copy of T»(£ — 1) for each variable node in the first 
generation and attach it to them. 

Assume that, for a typical random tree formula T(i), the marginal distribution of the variable at the 
root is weakly dependent on the values assigned at the boundary. Following the above lines, one can use 
this fact to prove that BP computes good approximations for the marginals in a random fc-SAT formula. 
In fact it turns out that an estimate of the formal 

E-r(t) sup \n(xi\x t ) - fJ,(xi\^)\ < e~ Kt (5.10) 

•£t '— t 

can be proved if the clause density a stays below a threshold a u (k) that is estimated to behave as 
a u (fc) = ii2sA[l + 0fc (l)]. 

While we refer to the original paper [3DJ for the details of the proof we limit ourselves to noticing 
that the left hand side of Eq. (|5.10[) can be estimated efficiently using a density evolution procedure. 
This allows to estimate the threshold a u (k) numerically. Consider in fact the log- likelihood (here we are 
identifying {True, False} with { + 1, —1}) 

2 M(-l|£t) 

This quantity depends on the assignment of the variables on the boundary, x t . Since we are interested 
on a uniform bound over the boundary, let us consider the extreme cases 

h (t) = max h (t \x t ) , h {t) = min h^(x t ) . (5.12) 

■Lt —t 

It is then possible to show that the couple (h^\h^) obeys a recursive distributional equation that, as 
mentioned, can be efficiently implemented numerically. 



6 Belief Propagation Beyond the Binary Symmetric Channel 

So far we have considered mainly the case of transmission over BMS channels, our reference example 
being the BSC. There are many other channel models that are important and are encountered in practical 
situations. Fortunately, it is relatively straightforward to extend the previous techniques and statements 
to a much larger class, and we review a few such instances in this section. 

6.1 Binary Memoryless Symmetric Channels 

In order to keep the notation simple, we assumed channel output to belong to a finite alphabet A. In our 

main example, the BSC, we had A = {0, 1}. But in fact all results are valid for a wider class of binary 

memoryless symmetric (BMS) channels. One can prove that there is no loss of generality in assuming 

the output alphabet to be the real line R (eventually completed with K = RU {±oo}). 

Let y = (yi, . . . , y^) be the vector of channel outputs on input x = (x\, . . . , xjf). For a BMS the 

input is binary, i.e. x £ {0, 1}^. Further the channel is memoryless, i.e. the probability density of 
— JV 

getting y S R at output when the input is x, is 

N 

Q(yk) = l[Q(yt\xt). 



Here the sup is taken over assignments x t , g/ t that can be extended to solutions of T(t). 
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Finally, the symmetry property can be written without loss of generality, as Q(yt\xt = 1) = Q{— yt\xt = 

0). ' ' ' 

One of the most important elements in this class is the additive white Gaussian noise (AWGN) 
channel, defined by 

y t = x t + z t , t E {1, . . . , N} , 

where the sequence {z t } is i.i.d. consisting of Gaussian random variables with mean zero and variance 
<T 2 . It is common in this setting to let xi take values in { + 1,-1} instead of {0, 1} as we have assumed 
so far. The AWGNC transition probability density function is therefore 

| (ij-x) 2 

Q(yt\x t ) 



V27R7 2 



The AWGNC is the basic model of transmission of an electrical signal over a cable (here the noise is 
due to thermal noise in the receiver) and it is also a good model of a wireless channel in free space (e.g., 
transmission from a satellite). 

Although the class of BMS channels is already fairly large, it is important in practice to go beyond 
it. The extension to the non-binary case is quite straightforward and so we will not discuss it in detail. 
The extension to channels with memory or the asymmetric case are more interesting and so we present 
them in the subsequent two sections. 



6.2 Channels With Memory 

Loosely speaking, in a memoryless channel the channel acts on each transmitted bit independently. In 
a channel with memory, on the other hand, the channel acts generally on the whole block of input bits 
together. An important special case of a channel with memory is if the channel can be modeled as a 
Markov chain, taking on a sequence of "channel states." Many physical channels posses this property 
and under this condition the message-passing approach can still be applied. For channels with memory 
there are two problems. First, we need to determine the capacity of the channel. Second, we need to 
devise efficient coding schemes that achieve rates close to this capacity. It turns out that both problems 
can be addressed in a fairly similar framework. Rather than discussing the general case we will look at 
a simple but typical example. 

Let us start by computing the information rate/capacity of channels with memory, assuming that 
the channel has a Markov structure. As we discussed in Section [2751 in the setting of BMS channels, the 
channel capacity can be expressed as the difference of two entropies, namely as H(X) — H(X\Y). Here, 
X denotes the binary input and Y denotes the observation at the output of the channel whose input is 
X. Given two random variables X, Y, this entropy difference is called the mutual information and is 
typically denoted by I(X; Y) = H(X) - H(X; Y). 

A general channel, is defined by a channel transition probability Q{y\ \xi) (here and below x^ 
denotes the vector (x\, . . . , xn)). In order to define a joint distribution of the input and output vectors, 
we have to prescribe a distribution on the channel input, call it p(x^). The channel capacity is obtained 
by maximizing the mutual information over all possible distributions of the input, and eventually taking 
the N — > oo limit. In formulae 

C(Q)= lim sup I(X^;Y 1 N )/N . 

For BMS channels it is possible to show that the maximum occurs for the uniform prior: p(x^) — 1/2 N . 
Under this distribution, I(X^ ; Y^)/N is easily seen not to depend on N and we recover the expression 
in Sec. [23] 

For channels with memory we have to maximize the mutual information over all possible distributions 
over {0, 1}^ (a space whose dimension is exponential in TV), and take the limit N — > oo. An easier task is 
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to choose a convenient input distribution p{ ■ ) and then compute the corresponding mutual information 
in the N — > oo limit: 

/= km HX^YfyN. (6.1) 

N— >oo 

Remarkably, this quantity has an important operational meaning. It is the largest rate at which we can 
transmit reliably across the channel using a coding scheme such that the resulting input distribution 
matches p( ■ ). 

To be definite, assume that the channel is defined by a state sequence {<7t}t>o, taking values in a 
finite alphabet, such that the joint probability distribution factors in the form 

n 

pW.I/i^o) =pM]Jp(xi,yi,a l | a^). (6.2) 

i=l 

We will further assume that the transmitted bits (x\, -,xn) are iid uniform in {0, 1}. The factor graph 
corresponding to (|6.2[) is shown in Fig. 1171 It is drawn in a somewhat different way compared to 
the factor graphs we have seen so far. Note that in the standard factor graph corresponding to this 
factorization all variable nodes have degree two. In such a case it is convenient not to draw the factor 
graph as a bipartite graph but as a standard graph in which the nodes correspond to the factor nodes 
and the edges correspond to the variable nodes (which have degree two and therefore connect exactly two 
factors). Such a graphical representation is also known as normal graph or as Forney-style factor graph 
(FSFG), in honor of Dave Forney who introduced them [15j . Let us now look at a concrete example. 



X2 £3 
(7 2 



0-3 


fn-1 






• • 1 


1 



p(a ) p(x2,y2,&2 I cri) 

p(xi,yi,cri I a Q ) p{x 3 , 7/3,0-3 I a 2 



p(x n ,y n ,(T n I 0"n-l) 

Figure 17: The FSFG corresponding to (|6,2p . 



Example 5: [Gilbert-Elliott Channel] The Gilbert-Elliot channel is a model for a fading channel, 
i.e., a channel where the quality of the channel is varying over time. In this model we assume that 
the channel quality is evolving according to a Markov chain. In the simplest case there are exactly 
two states, and this is the original Gilbert-Elliott channel (GEC) model. More precisely, consider 
the two-state Markov chain depicted in Fig. [T51 Assume that {X t } t >i is i.i.d., taking values in 
{±1} with uniform probability. 




Figure 18: The Gilbert-Elliott channel with two states. 



The channel is either in a good state, denote it by G, or in a bad state, call it B. In either state 
the channel is a BSC. Let the crossover probability in the good state be cq and in the bad state 
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be £b, with < ec < £b < 1/2. Let P be the 2x2 matrix 

"-{',1 

which encodes the transition probabilities between the states (the columns indicate the present 
state and the rows the next state). Define the steady state probability vector p = (pg,Pb), i.e., the 
vector which fulfills Pp T — p T . This means that in steady state the system spends a fraction pa of 
the time in state G and a fraction pe of the time in state B. If we consider e.g. the state G then 
the detailed balance condition reads peg = p-ob. From this we get p — (b/ (g + b), g/(g + b)). More 
generally, let us assume that we have s states, s € N, and that the channel in state i, i £ [s], is 
the BSC(ei). Let P be the s x s matrix encoding the transition probabilities between these states. 
Let p denote the steady-state probability distribution vector. If (7 — P T + E) is invertible then a 
direct check shows that p — e(I — P T + E) -1 , where e is the all-one vector of length s, I is the 
s x s identity matrix and E is the s x s all-one matrix. 

Note that the state sequence is ergodic as long as the Markov chain is irreducible (i.e. there is a 
path of strictly positive probability from any state to any other state) and aperiodic (i.e. there 
there exists such a path for any number of steps large enough). In the original Gilbert-Elliot model 
this is true as long as < g, b < 1. 

Consider the computation of the maximal rate at which we can transmit reliably. We have 

I(XF;Y?) = H{Y?) - H{Y?\X™). 

Let us see how we can compute limAr^oo H(Y^)/N . Because of the ergodicity assumption on the state 
sequence, —4- logp(yj v ) converges with probability one to limjv^oo H(Y-j N )/N. It follows that if we can 
compute — jy logp(yj v ) for a very large sequence, then with high probability the value will be close to 
the desired entropy rate. Instead of computing p{y^), let us compute p((Jn, Vi)- From this we trivially 
get our desired quantity by summing, 

p(vf) = X>(<™,yi w ). 

Note that 

P{vN,Vl)= ^2 P(xN,<TN-l,VN,yi) 
KNi&N—l 

= T~] p(x N ,a N ,y N \a N ^i)p{(j N ^ 1 ,y^~ 1 ) . (6.3) 

v ' S v ' 

kernel message 

From this we see that p{<jniUi) can De computed recursively. In fact this recursion corresponds to 
running the BP message-passing rules on the factor graph depicted in Fig. 1171 (which is a tree): denote 
the message which is passed along the edge labeled by crjv by vn{pn)- Then according to the BP 
message-passing rules we have 

vn(&n) — p(xn, <tn, Vn I Cjv-l) Vn-x{(TN-\)- 

If we compare this to the recursion stated in (|6.3p we see that these two recursions are identical. In 
other words, vn{vn) = p(&n, V\), so that 

lim H{Y 1 N )/N = - lim log(Vi/jv(crjv)W. (6-4) 
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From a practical perspective it is typically more convenient to pass normalized messages vn(o~n) so that 
vn(o~n) = 1- The first message Vq((Tq) — p(o~o) is already a probability distribution and, hence, 
normalized, vo(ao) — ^o(cxo)- Compute vi(ai) and let Ai = vi(<j\). Define v\(ai) = vi{<t\)/\\. 
Now note that by definition of the message-passing rules all subsequent messages in the case of rescaling 
differ from the messages which are sent in the unsealed case only by this scale factor. Therefore, if 
denotes the normalization constant by which we have to divide at step i so as to normalize the message 
then Dn{&n) — ^n^n) / (YiiLi ^i)- It follows that 

lim H(Y 1 N )/N = - lim logfV a N (a N ))/N 

& N 

N N 

It remains to compute H(Yf | ). We write H(Yf \ X?)/N = H(Y 1 N ,X^)/N - H(X^)/N. 
The second part is trivial since the inputs are i.i.d. by assumption so that H(X^)/N = 1. For the 
term H{Y^ ,X^)/N we use the same technique as for the computation of H(Y^)/N. Because of 
the ergodicity assumption on the state sequence, — logp(yf r , xf) converges with probability one to 
limjv^oo H(Y^ , Xi)/N . We write p{yi — J2a N p(°~n, Ui , Xi) and use the factorization 

p((T N ,y 1 ,X{) = P[CTN-l,CTN,Vl ,Xl ) 

= y~] p{x N ,o- N ,y N | ctat-i) ■ p(o- N -i,yi ~ X ,Xi _1 ) . 

V v ' V v ' 

kernel message 

In words, we generate a random instance Xf and and run the BP algorithm on the FSFG shown 
in Fig. [T7| assuming that both Yf and Xf are 'quenched.' Taking the logarithm, multiplying by minus 
one and normalizing by 1 /N gives us an estimate of the desired entropy. 

Now that we can compute the maximal rate at which we can transmit reliably, let us consider coding. 
The symbol MAP decoder is 

Xi(y) = argma,x Xi p(xi | ) 

El N N N\ 
P\ x i ,Vi ,<?o ) 

N 

= argmax^ p{a Q ) p{xj , Vj, <Jj I CTj_i)I £ (xf ). 

{x S ,jjti} j=l 

In words, the FSFG in Fig. \T7\ describes also the factorization for the message-passing decoder if we add 
to it the factor nodes describing the definition of the code. As always, this factor graph together with 
the initial messages stemming from the channel completely specify the message-passing rules, except for 
the message-passing schedule. Let us agree that we alternate one round of decoding with one round of 
channel estimation. No claim as to the optimality of this scheduling rule is made. 

Notice that the correlations induced by the markovian structure of the channel arc in general short 
ranged in time. This is analogous to what happens with a one-dimensional spin model, whose correlation 
length is always finite (at non-zero temperature). A good approximation to the above message passing 
schedule is therefore obtained by a 'windowed' decodes. This means that the state at time t is estimated 
only of the basis of observations between time t — R and t + R, for some finite R. 

Assuming windowed decoding for channel estimation, it is not hard to show that after a fixed number 
of iterations, the decoding neighborhood is again asymptotically tree-like. In the case of the GEC the 
channel symmetry can be used to reduce to the all-zero codeword. Therefore, we can employ the 
technique of density evolution to determine thresholds and to optimize the ensembles. 



36 



































1 






















































































































































































I 



























-10 10 20 30 40 50 60 70 



0.3- 
0.2 
0. 







































































































































1 






























































































) 































-10 10 20 30 40 50 60 70 

















3- 














2- 


±: 












1- 















-10 10 20 30 40 50 60 70 



0.3 
0.2 
0. 



-10 10 20 30 40 50 60 70 



























































































































5-4-3-2-1 1 i 


I : 




L 5 



























































































































5-4-3- 


2-1 o : 


L : 


) ; 


5 - 


L 5 

































































































I 










A . . 










1 


5-4-3-2-1 : 


l 2 ; 


3 - 


I 5 







































































































i ■ 






L 







-5-4-3-2-1 12 3 4 5 



Figure 19: Density evolution for the GEC at iteration 1, 2, 4, and 10. The left pictures show the densities 
of the messages which are passed from the code towards the part of the FSFG which estimates the channel 
state. The right hand side shows the density of the messages which are the estimates of the channel state 
and which are passed to the FSFG corresponding to the code. 



Example 6: [GEC: State Estimation] For the case of transmission over the GEC the iterative 
decoder implicitly also estimates the state of the channel. Let us demonstrate this by means of 
the following example. We pick a GEC with three states. Let 

/ 0.99 0.005 0.02 \ 
P = 0.005 0.99 0.02 , 
\ 0.005 0.005 0.96 / 

which has a steady state probability vector of p sa (0.4444, 0.4444, 0.1112). Finally, let the channel 
parameters of the BSC in these three states be (61,62,63) « (0.01,0.11,0.5). This corresponds to 
an average error probability of e avg = Y^i=\Pi e i ~ 0-108889. Using the methods described above, 
the capacity of this channel (assuming uniform inputs) can be computed to be equal to C ~ 0.583 
bits per channel use. This is markedly higher than 1 — h(e avg ) « 0.503444, which is the capacity 
of the BSC(e avg ). The last channel is the channel which we experience if we ignore the Markov 
structure. 

Fig. [T5] shows the evolution of the densities for an optimized ensemble of rate r « 0.5498. The 
pictures on the right correspond to the messages which are passed from the part of the factor 
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graph which estimates the state towards the part of the factor graph which describes the code. 
These messages therefore can be interpreted as the current estimate of the state the channel is 
in at a given point in time. Note that after 10 iterations 5 clear peaks emerge. These peaks are 
at ±log(0.99/0.01) fa ±4.595, ± log(0.9/0.1) « ±2.197, ± log(0.5/0.5) = 0. They correspond to 
the received likelihoods in the three possible channel states. In other words, the emergence of the 
peaks shows that at this stage the system has identified the channel states with high reliability. 
This is quite pleasing. Although the channel state is not known to the receiver and can not be 
observed directly, in the region where the iterative decoder works reliably it also automatically 
estimates the channel state with high confidence. 

Although we only looked a very particular example it is quite typical of the general situation: as 
long as the channel memory can be described by a Markov chain the factor graph approach applies and 
we can use message-passing schemes to construct efficient coding schemes [TH HI [211 E21 EH [43] . 

6.3 Asymmetric Channels - The Z Channel 

Let us now consider the second generalization, namely the case of non- symmetric channels. 

Consider the channel depicted on the right of Fig. [TJ For obvious reasons it is called the Z channel 
(ZC). This channel has binary input and it is memoryless but it is not symmetric. Nevertheless, essen- 
tially the same type of analysis which we performed in Section [4] can be applied to this case as well. 
Symmetry is therefore a nice property to have but it is not essential. 

Consider the capacity of this channel. Since the channel is not symmetric the capacity is not neces- 
sarily given by the mutual information between channel input and channel output for a uniform input 
distribution of the input. We must instead maximize the mutual information over the input distribu- 
tion. Since the channel is memoryless, it can be assumed that the input is given by a sequence of i.i.d. 
bernoulli variables. Assuming that p{xi = 0) = a, the output distribution is 

(p(Vi = 0),p{yi = 1)) = (ap, 1 - ap), 

so that the mutual information I a (X;Y) for a fixed a is equal to 

I a (X; Y) = H(Y) - H(Y \ X) = h(ap) - ah(p). (6.5) 

Some calculus reveals that the optimal choice of a is 

pP/P 

^ ' 1 + ppP/P 

so that 

Czc(p) = h(a(p)p) - a(p)h(p). 

Fig. [20] compares Czc(p) with I a _i(X;Y), i.e., it compares the capacity with the transmission rate 
which is achievable with uniform input distribution. This is important and surprising - only little is 
lost by insisting on an uniform input distribution: the rate which is achievable by using a uniform input 
distribution is at least a fraction |eln(2) ps 0.924 of capacity over the entire range of p (with equality 
when p approaches one) . Even more fortunate, from this perspective the Z channel is the extremal case 
[321 154] : the information rate of any binary-input memoryless channel when the input distribution is 
the uniform one is at least a fraction ieln(2) of its capacity. From the above discussion we conclude 
that, when dealing with asymmetric channels, not much is lost if we use a binary linear coding scheme 
(inducing a uniform input distribution). 

Consider the density evolution analysis. Because of the lack of symmetry we can no longer make 
the all-one codeword assumption. Therefore, it seems at first that we have to analyze the behavior of 
the decoder with respect to each codeword. Fortunately this is not necessary. First note that, since we 
consider an ensemble average, only the type of the codeword matters. More precisely, let us say that a 
codeword has type r if the fraction of zeros is r. For x € £, let t(x) be its type. Let us assume that we 
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Figure 20: Comparison of Czc(p) (solid curve) with I a _i(X;Y) (dashed curve), both measured in bits. 



use an LDPC ensemble whose dominant type is one-half. This means that "most" codewords contain 
roughly as many zeros as one. Although it is possible to construct degree-distributions which violate 
this constraint, "most" degree distributions do fulfill it. Under this assumption there exists some strictly 
positive constant 7 such that 

P {r{x) £ [I/ 2 - S/y/n, 1/2 + < e~ 5 ^ , (6.7) 

where the probability is with respect to a uniformly random codeword x. We can therefore analyze the 
performance of such a system in the following way: determine the error probability assuming that the 
type of the transmitted codeword is "close" to the typical one. Since sublinear changes in the type do 
not figure in the density analysis, this task can be accomplished by a straightforward density evolution 
analysis. Now add to this the probability that the type of a random codeword deviates significantly 
from the typical one. The second term can be made arbitrarily small (see right hand side of (|6.7j) ) by 
choosing 8 sufficiently large. 

We summarize: if we encounter a non-symmetric channel and we are willing to sacrifice a small 
fraction of capacity then we can still use standard LDPC ensembles (which impose a uniform input 
distribution) to transmit at low complexity. If it is crucial that we approach capacity even closer, a more 
sophisticated approach is required. We can combine LDPC ensembles with non-linear mappers which 
map the uniform input distribution imposed by linear codes into a non-uniform input distribution at 
the channel input in order to bring the mutual information closer to capacity. For a detailed discussion 
on coding for the Z-channel we refer the reader to [34] [59l [5] . 

7 Open Problems 

Let us close by reviewing some of the most important open challenges in the channel coding problem. 

7.1 Order of Limits 

Density evolution computes the limit 

lim lim EpP^l. 

In words we determined the limiting performance of an ensemble under a fixed number of iterations as 
the blocklength tends to infinity and then let the number of iterations tend to infinity. As we have seen, 
this limit is relatively easy to compute. What happens if the order of limits is exchanged, i.e., how does 
the limit 

lim lim EpP'* 1 ] 
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behave? This limit is closer in spirit to the typical operation in practice: for each fixed length the BP 
decoder continues until no further progress is achieved. We are interested in the limiting performance 
as the blocklength tends to infinity. 

For the BEC it is known that the two limits coincide. If we combine this with the fact that for the 
BEC the performance is a monotone function in the number of iterations (any further iteration can only 
make the result better) then we get the important observation that regardless of how we take the limit 
(jointly or sequentially), as long as both the blocklength as well as the number of iterations tend to 
infinity we get the same result. From a practical perspective this is comforting to know: it shows that 
we can expect a certain robustness of the performance with respect to the details of the operation. 

It is conjectured that the same statement holds for general BMS channels. Unfortunately, no proof 
is known. 

7.2 Finite-Length Performance 

The threshold gives an indication of the asymptotic performance: for channel parameters which are 
better than the threshold sufficiently long codes allow transmission at an arbitrarily low probability of 
bit error. If, on the other hand, we transmit over a channel which has a parameter that is worse than the 
threshold then we can not hope to achieve a low probability of error. This is an important insight but 
from a practical perspective we would like to know how fast the finite-length performance approaches this 
asymptotic limit. There can be many different ensembles that all have the same asymptotic performance 
but that might have a substantially different finite-length behavior. Can we predict which one we should 
choose a priori without having to resort to simulations? The typical convergence of the performance to 
the asymptotic limit is shown in Fig. 1211 The points correspond to simulation results whereas the solid 
curves correspond to a general scaling conjecture [2]. Roughly speaking, this scaling conjecture states 
that around the threshold the error probability behaves as follows: Let the channel be parameterized 
by e with increasing e indicating a worsening of the channel. Let £d be the BP threshold, and define 
z = \/n(e — ed). Then for z fixed and n tending to infinity we have 

P B (n,e) = $(z/a) (l + o(l)), 

where $( • ) is the error function (i.e. $(x) is the probability that a standard normal random variable is 
smaller than x), and a is a constant which depends on the channel as well as on the channel. 

For the BEC this scaling law has been shown to be correct [3J. In fact, even a refined version is 
known [T3]: define z = y/n(e — ed + (3n~^) where (3 is a constant depending on the ensemble. Then for 
z fixed and n tending to infinity we have 

P B (n,e) = Q(z/a)(l + 0(n- 1 / 3 )). 

For general channels on the other hand the problem is largely open. If proved to be correct, finite 
length scaling laws could be used as a tool for an efficient finite-length optimization. 

7.3 Capacity- Achieving Codes 

For the most part we have taken the point of view that we are given an ensemble of codes and a family 
of channels and we would like to determine the performance of this combination. For instance, the 
most fundamental question is to determine the threshold noise for such a code/decoder combination. 
Hopefully this threshold is close to the best possible as determined by Shannon's capacity formula. 

But we can take a more active point of view. Given a family of channels how should we choose the 
ensemble in order for the threshold noise to be as high as possible. In other words, can we approach the 
capacity of the channel? 

For the BSC this question has been answered by Luby, Mitzenmacher, Shokrollahi, Spielman, and 
Steman in [27]. These authors showed that by a suitable choice of the degree distribution one can 
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Figure 21: Performance of the BP decoder for the (3, 6) -regular ensemble when transmission takes place 
over the BSC. The block lengths are n = 2 l , i = 10, • • ■ , 20. The dots correspond to simulations. For most 
simulation points the 95% confidence intervals are smaller than the dot size. The lines correspond to the 
analytic approximation based on scaling laws. 



approach the capacity arbitrarily closely. More precisely, in order to approach capacity up to a fraction 
5 the average degree has to grow like log(l/<5) and this is the best possible. For general channels it is 
not known whether capacity can be achieved. Although the resolution of this problem will most likely 
only have a small practical implication it is without doubt the most important open theoretical problem 
in this area. 

A A Generating Function Calculation 

We want to compute the number of ways of picking n distinct objects from M groups each containing 
k objects, in such a way that the number of selected elements in each group is even. Let us denote this 
number by C n . The result of this computation is used in Section [3. 3[ where we claimed the result to be 
C n = coeff[qk{z) M , z n ] (in that case n = Iw). 

First notice that given mi, ... ,mM all even, the number of ways of choosing mi objects from the 
first group, mi from the second, etc is 




(A.l) 



The desired number C n is the sum of this quantity over the even numbers mi,-- - ,m n such that 
mi, . . . , mM G {0, . . . k} and mi + • • • + Tom = n. The corresponding generating function is 

n— mi...m n v ' ' 

the sum being restricted over even integers mi, . . . ,mu € {0, ...k}. We now notice that the sum 
factorizes yielding 




(A.3) 



which proves the claim. 
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